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SECTION  1 
INTRODUCTION 


This  report  documents  the  numerical  algorithms  in  the  EPIC  hydrocode.  The  EPIC  code  was 
conceived  in  the  early  1970s  and  has  imdergone  significant  development  since  that  time.  The 
first  two  technical  publications  for  the  2D  and  3D  versions  of  EPIC  are  given  in  References  1  and 
2,  and  the  first  two  contract  reports  are  given  in  References  3  and  4.  Although  most  of  the 
numerical  algorithms  in  EPIC  have  been  published  as  technical  papers  and/or  reports,  there  has 
not  been  a  recent  publication  that  contains  all  of  the  important  algorithms  within  a  single  report. 

In  Section  2  of  this  report  there  is  a  description  of  the  structure  of  the  EPIC  code.  This  is 
included  to  show  where  and  how  the  numerical  algorithms  fit  into  the  computational  framework. 

Section  3  represents  the  main  portion  of  this  report  and  contains  all  of  the  finite  element 
algorithms.  This  contains  several  ID,  2D  and  3D  geometry  options,  and  each  of  these  geometries 
contains  multiple  element  types.  The  contact,  sliding  and  erosion  algorithms  are  also  included  in 
this  section. 

The  SPH  (Smooth  Particle  Hydrodynamics)  algorithms  are  presented  in  Section  4,  and  the 
linking  together  of  SPH  nodes  and  standard  finite  elements  is  presented  in  Section  5. 

The  report  concludes  with  Section  6,  which  describes  the  numerous  material  models  available  in 
EPIC. 
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SECTION  2 

EPIC  CODE  STRUCTURE 


This  section  provides  an  overview  of  the  structure  of  the  EPIC  code.  Figure  1  shows  a  hierarchy 
chart  of  some  of  the  key  subroutines  in  the  EPIC  code.  Because  there  are  hundreds  of 
subroutines  in  the  code  it  is  not  possible  to  include  all  of  them  in  a  single  figure.  In  general, 
however,  the  flow  of  the  code  can  be  described  by  Figure  1.  In  most  instances  the  calling 
sequence  goes  from  left  to  right. 

The  following  provides  a  brief  description  of  the  fimctions  of  the  subroutines  shown  in  Figure  1 : 
EPIC  Calls  the  first  layer  of  subroutines. 

HD  AT  A  Reads  the  input  file. 

GEOM  Generates  the  initial  geometry. 

MASS  Computes  nodal  masses. 

NVECT  Groups  node  arrays  into  blocks  to  allow  for  vectorized  computational  loops. 

START  Assigns  initial  velocities  and  explosive  detonation  points  and  times. 

HEAT  Computes  nodal  heat  conduction  quantities  from  element  heat  conduction 

parameters. 

WRITEG  Writes  output  file  for  the  Preprocessor  and  computes  initial  integration  time 

increment. 

CHUNK  Computes  data  (velocity,  position,  energy,  momentum,  etc.)  for  user- 

specified  chunks  of  elements. 

SDATA  Computes  data  (velocity,  position,  energy,  momentum,  etc.)  for  the  entire 

system. 
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Figure  1.  Partial  Hierarchy  Chart  for  EPIC  Subroutines 


RECALL 


Reads  the  restart  file  for  restart  runs. 


LOOP 

SAVE 

NLOOP 

MOTION 

MOVE 

SLIDEl 

SLIDE2 

SLIDES 

ELOOP 

EGET 

VOLUME 

GMCON 


Calls  the  three  primary  subroutines  for  each  computational  cycle. 

Writes  the  restart  file  for  subsequent  postprocessing  and  restart  runs. 

Calls  the  subroutines  that  govern  the  velocities  and  positions  of  the  nodes. 
This  is  the  first  of  three  primary  subroutines/fimctions  in  the  computational 
cycle. 

Updates  the  velocities  of  the  nodes  through  the  equations  of  motion. 

Updates  the  displacements  of  the  nodes  through  the  equations  of  motion. 

Calls  subroutines  to  update  the  velocities  and  displacements  of  the  nodes  for 
ID  contact  interfaces. 

Calls  subroutines  to  update  the  velocities  and  displacements  of  the  nodes  for 
2D  sliding  interfaces. 

Calls  the  subroutines  to  update  the  velocities  and  displacements  of  the  nodes 
for  3D  sliding  interfaces. 

Calls  the  subroutines  that  compute  the  element  quantities.  This  is  achieved 
by  computing  with  blocks  of  elements  that  have  the  same  material  and 
element  type.  This  is  the  second  of  the  primary  subroutines/fimctions  in  the 
computational  cycle. 

Gathers  nodal  data  (velocities  and  displacements)  and  transforms  it  into 
element  data  such  that  the  element  computations  can  be  vectorized. 

Computes  volumetric  strains  and  strain  rates. 

Computes  geometry  factors  that  are  used  for  strain  rate,  force  and  heat 
conduction  computations. 
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STRAIN 


Computes  normal  and  shear  strain  rates. 


STRESS 

SOLID 

SHELLS 

MEGRU 

ARTVIS 

HEBURN 

CRUSH 

BRITLE 

REACT 

DAMAGE 

FAIL 

CORANT 


Calls  the  subroutines  that  compute  the  stresses  in  the  elements. 

Computes  the  maximum  allowable  Von  Mises  stress  for  most  of  the 
material  models  that  exhibit  strength.  Also  calls  some  additional 
subroutines  for  more  complicated  strength  models. 

Performs  iterative  strain  rate  and  stress  computations  for  shells  with  bending 
capability. 

Computes  pressure,  internal  energy  and  sound  velocity  for  solids  (metals). 
Computes  artificial  viscosity. 

Computes  pressure,  internal  energy  and  soimd  velocity  for  explosives  that 
use  a  programmed  bum  algorithm. 

Computes  pressmre,  internal  energy  and  sound  velocity  for  crushable 
(concrete)  materials. 

Computes  pressure,  internal  energy  and  sound  velocity  for  brittle  (ceramic) 
materials. 

Calls  subroutines  to  compute  bum  fractions,  pressure  and  sound  velocity  for 
reactive  explosives. 

Computes  damage  for  fracture. 

Totally  fails  highly  distorted  elements  that  have  exceeded  a  user-specified 
strain. 

Computes  the  allowable  integration  time  increment  for  the  next 
computational  cycle. 
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FORCE 

CDUCTl 

EPUT 

CDUCT2 

NALOOP 


Computes  nodal  forces  that  result  from  the  stresses  in  the  elements. 

Computes  heat  flow  across  the  elements  for  the  heat  conduction  option. 

Scatters  the  element  forces  to  the  corresponding  nodes. 

Redefines  element  temperatures  and  internal  energies  based  on  temperatures 
and  internal  energies  of  the  nodes.  Used  only  for  the  heat  conduction 
option. 

Calls  the  subroutines  that  compute  SPH  nodal  quantities.  This  is  the  third  of 
the  primary  subroutines/fimctions  in  the  computational  cycle.  The 
subroutines  called  by  NALOOP  are  similar  to  those  called  by  ELOOP. 
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SECTION  3 

FINITE  ELEMENT  ALGORITHMS 


This  section  provides  a  description  of  the  algorithms  for  the  finite  elements  in  the  EPIC  Code. 
For  ID  geometry,  the  algorithm  for  Cartesian  geometry  is  presented  in  detail.  The  other  ID 
options  for  cylindrical  and  spherical  geometry  only  include  descriptions  of  those  features  that  are 
different  from  the  Cartesian  geometry. 

The  baseline  2D  geometry  is  the  axisymmetric  case,  and  this  is  presented  in  detail.  The  other  2D 
geometries  (axisymmetric  with  spin,  plane  strain  and  plane  stress)  only  include  descriptions  of 
those  features  that  are  different  from  the  axisymmetric  geometry.  There  is  only  one  geometric 
option  for  3D  geometry  and  it  is  presented  in  detail.  The  material  models  to  determine  stresses 
are  not  presented  in  this  section,  but  are  presented  in  Section  6. 

An  attempt  has  been  made  to  keep  the  nomenclature  consistent  for  the  various  algorithms,  but 
this  is  not  always  possible  because  this  work  has  been  performed  by  several  different  people  over 
the  course  of  many  years.  The  general  philosophy  is  to  have  the  nomenclature  in  this  report  be 
similar  to  that  which  exists  in  the  code. 

A  final  comment  is  that  the  numbering  of  the  equations  begins  with  Equation  1  at  the  beginning 
of  each  subsection. 

3.1  ID  CARTESIAN  GEOMETRY 

A  description  of  the  ID  Cartesian  geometry  is  shown  in  Figure  2.  The  z  coordinate  is  the 
coordinate  of  interest  and  the  element  is  defined  by  nodes  i  and  j,  with  node  i  having  a  higher  z 
coordinate  than  node  j  (zi  >  zj).  The  algorithm  is  for  ID  imiaxial  strain  (not  uniaxial  stress); 
therefore  the  cross-sectional  area  is  always  A  =  1.0. 

The  mass  at  node  i,  for  an  individual  element,  is 

M-PoV„/2  .1. 
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ID  Cartesian  Element 


Assemblage  of  Elements 


Contact 

Conditions 


Z 
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Figure  2.  Description  of  the  1D  Cartesian  Element 
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where  po  is  the  initial  density  of  the  material  and  the  initial  volume  is 


Vo  =  A(Zi-z.)  =  Zi-z. 

forA=  1.0. 


(2) 


When  element  C  is  incorporated  into  an  assemblage  of  elements,  as  shown  in  Figure  2,  then  the 
total  mass  at  node  i  contains  the  individual  element  masses  from  both  elements  B  and  C.  The 
total  mass  at  node  i  is 

(3) 


3.1.1  Equations  of  Motion 


The  acceleration,  velocity  and  position  of  node  i  are  determined  as  follows: 
z|  =  Fj/Mj 

zr=zr+zi^ 

=z|+z‘^At 


(1) 

(2) 

(3) 


In  Equation  1  the  acceleration  is  z|  at  time  =  t.  F^'  is  the  net  z  direction  force  on  node  i,  and  this 
is  the  sum  of  the  individual  element  forces  on  node  i.  These  forces  are  obtained  from  the 
element  computations  in  the  previous  cycle  of  integration.  This  will  be  explained  in  greater 
detail  in  subsection  3.1.3.  Mj  is  the  total  mass  at  node  i,  as  given  in  Equation  3  of  subsection 

3.1. 


The  updated  velocity  in  Equation  2  is  constant  for  the  interval  between  time  =  t  and  time  =  t  +  At. 
The  constant  velocity  for  the  previous  time  increment  is  z\~  and  At  is  the  average  of  the  two 

integration  time  increments  about  time  =  t. 


The  integration  time  increment  is  limited  to 
At  =  C,[h/(^  +  Vg^+c^ ) 


(4) 
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where  =  CqQ/p,  h  =  Zi  -  zj  is  the  height  of  the  element  and  Cs  is  the  sound  velocity  of  the 
material  (References  1  and  5).  Q  is  the  artificial  viscosity  and  Cq  is  a  constant  for  the  artificial 
viscosity.  These  are  described  in  Section  6.  The  Courant  sound  speed  fraction,  Q,  must  be  less 
than  unity  (Q  <  1.0)  to  ensure  numerical  stability.  Q  =  0.9  is  a  typical  value. 

3.1.2  Contact  Interfaces 

Contact  interfaces  can  be  handled  in  a  straightforward  manner  in  ID.  The  user  must  specify  two 
nodes  for  each  interface.  Node  s  is  designated  as  the  slave  node  and  it  is  the  node  above  (higher 
z  coordinate)  master  node  m.  The  slave  and  master  designations  are  used  for  convenience  and  do 
not  cause  a  bias  in  the  results. 

After  the  equations  of  motion  have  been  applied,  a  check  is  made  to  determine  if  there  is 
interference  (zm  ^  Zg).  If  there  is  no  interference,  then  there  is  no  contact  and  the  surfaces  are 

free. 


If  there  is  contact,  it  is  necessary  to  adjust  the  positions  and  velocities  of  the  two  nodes.  The 
center  of  mass  of  the  two  nodes  is 


z  =  (m,z,  +  M„,z„  )  /  (m,  +  M„  ) 

(1) 

where  Mg  and  Mm  are  the  masses  of  the  interface  nodes,  and  Zg 

coordinates. 

and  Zm  are  the  corresponding  z 

The  movements  of  the  two  nodes,  to  the  center  of  mass  are 

Az3  =  z-z, 

(2) 

Az„  =  z-z„ 

(3) 

and  the  corresponding  velocity  changes  are 

AZj  =  AZj  /  At 

(4) 

6 

II 

5 

(5) 
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The  updated  positions  and  velocities  then  become 


„new 

Zs 


=  zf +AZ3 


(6) 


^new 

Zm 


,old 


+  Az„ 


(7) 


•  new 

Zs 


+  Az, 


(8) 


•new 


+  ^m 


(9) 


This  formulation  conserves  momentum  and  provides  consistency  between  the  velocities  and 
positions.  If  there  was  contact  during  the  previous  cycle,  then  the  velocities  of  the  two  contact 
nodes  are  identical.  If  there  was  a  gap  during  the  previous  cycle  then  there  will  continue  to  be  a 
closing  velocity  (but  only  for  one  cycle). 

3.1.3  ID  Cartesian  Elements 


There  are  three  primary  features  to  all  of  the  element  algorithms.  The  first  is  to  obtain  the  strains, 
strain  rates  and  rotational  rates.  From  these,  the  stresses  and  pressures  are  obtained  through  a 
variety  of  material  models.  Then,  equivalent  nodal  forces  must  be  determined  from  the  stresses. 

The  volumetric  strain  and  the  volumetric  strain  rate  are  as  follows: 

e,=V/V„-l  (1) 

e,=(er-<)/At  (2) 

In  Equation  1,  V  and  V©  are  the  current  and  initial  element  volumes;  and  in  Equation  2,  and 
are  the  volumetric  strains  at  time  =  t  +  At  and  time  =  t.  The  volumetric  strains  and  strain  rates 
are  based  on  the  initial  configuration,  as  opposed  to  the  shear  and  deviator  strain  rates  which  are 
based  on  the  current  configuration. 
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A  more  common  notation  is  to  represent  the  z  velocity  by  v ,  instead  of  z .  The  z  velocities  in 
the  element  are  then  assumed  to  vary  in  a  linear  manner,  such  that 

v  =  aj+a2z  (3) 

where  ai  and  a2  are  geometry  and  velocity  constants. 

Substituting  the  two  nodal  velocities  (vj  and  Vj)  and  the  two  nodal  positions  (zf  and  zj)  into 

Equation  3  gives  two  equations  and  two  unknowns  (ai  and  02).  The  strain  rate  in  the  z  direction 
can  then  be  determined  from 

All  of  the  other  strain  rate  components  are  zero. 

ex=ey=Yxy=Yxz=tyz=o  (5) 

where  and  Ey  are  the  other  two  normal  strain  rates  and  Y^^y ,  Yxz  Yyz  shear 

strain  rates. 

Some  of  the  constitutive  models  require  an  equivalent  strain  rate,  and  this  is  expressed  as 


which  reduces  to  e  =  /  3  for  the  special  case  of  ID  Cartesian  geometry. 

The  constitutive  models  also  require  deviator  strain  rates,  which  are  expressed  as 


^ave 

(7) 

(8) 

^ave 

(9) 
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where  Sx,  Sy  and  Sz  are  the  normal  deviator  stresses,  P  is  the  hydrostatic  pressure  and  Q  is  the 
artificial  viscosity.  These  are  described  in  detail  in  Section  6.  The  shear  stresses  are  txy  =  Txz  = 
Tyz  =  0  for  ID  geometry. 

Trial  values  of  the  deviator  stresses  at  time  =  t  +  At  are 


=  s^  +2Ge^At 

(13) 

=  Sy +2GeyAt 

(14) 

t+At 

“z 


=  S2  +  2Ge^At 


(15) 


In  Equation  13,  the  first  term  (s‘^)  is  the  deviator  stress  at  the  previous  time  and  the  second  term 
(2Ge,^At)  is  the  incremental  stress  due  to  the  incremental  strain  (e^^At)  during  that  time 
increment,  where  G  is  the  elastic  shear  modulus. 


Equations  13-15  assume  an  elastic  response  of  the  material.  If  the  strength  of  the  material  is 
exceeded,  then  plastic  flow  (or  fracture)  will  occur.  The  Von  Mises  yield  criterion  is  used  to 
determine  an  equivalent  stress,  a ,  that  can  be  compared  to  the  uniaxial  tensile  (or  compressive) 
strength  of  the  material.  The  general  form  of  the  equivalent  stress  is 
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Using  deviator  stresses  (instead  of  total  stresses)  and  setting  the  shear  stresses  to  zero  (for  ID 
response),  the  equivalent  stress  can  be  simplified  to 


If  a  is  not  greater  than  the  equivalent  tensile  strength  of  the  material,  a,  the  final  deviator  and 
shear  stresses  are  as  given  in  Equations  13—15.  If  o  is  greater  than  cr,  then  the  stresses  in 
Equations  13-15  are  multiplied  by  the  factor  (a  /  ct)  .  When  the  reduced  deviator  and  shear 
stresses  are  put  into  Equation  17,  the  result  is  always  a  =  a .  This  is  known  as  the  radial  return 
algorithm.  The  various  material  strength  models  for  c  are  presented  in  Section  6. 

During  plastic  flow  it  is  sometimes  necessary  to  determine  the  equivalent  plastic  strain  for  strain 
hardening  effects  on  the  strength  of  the  material,  or  to  determine  if  the  material  has  failed.  The 
first  step  in  this  process  is  to  adjust  the  total  strain  rates  to  plastic  strain  rates  by  subtracting  out 
the  elastic  portion  of  the  strain  rates. 


(18) 

e;  =  e,-(s“-s;)/2GAt 

(19) 

e;=e.-(sr-s:)/2GA« 

(20) 

Again,  the  plastic  shear  strain  rates  are  zero  for  ID  geometry.  The  expression  for  the  ID 
equivalent  plastic  strain  rate  is 


(21) 


The  equivalent  plastic  strain,  Ep ,  is  then  obtained  by  integrating  with  respect  to  time: 


e‘-"^‘=e;+EpAt 

After  the  element  stresses  are  obtained,  it  is  necessary  to  determine  concentrated  forces  to  act  on 
the  concentrated  masses  at  the  nodes.  These  nodal  forces  are  used  in  the  equations  of  motion  for 
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the  subsequent  integration  cycle.  The  net  z  forces  on  nodes  i  and  j,  for  an  individual  element,  are 
simply 


f;=-(5,a=-o. 

(23) 

F'  =  o,A  =  a. 

(24) 

for  A  =  1.0 .  Note  that  =  0  for  each  element,  and  this  ensures  equilibrium  for  the  system. 

This  results  in  conservation  of  momentum  if  there  are  no  external  forces  or  restraints.  The  final 
net  force  on  node  i  is 

(25) 


It  is  also  possible  to  add  external  forces  through  applied  pressures. 

The  heat  conduction  algorithm  assumes  the  temperature  varies  linearly  between  nodes  such  that 

T  =  03+047  (26) 

where  T  is  the  temperature  in  the  element,  and  O3  and  O4  are  geometry  and  nodal  temperature 
constants. 

Substituting  the  two  nodal  temperatures  (Tj  and  Tj)  and  the  two  nodal  positions  (zj  and  zj)  into 
Equation  26  gives  two  equations  and  two  unknowns  (O3  and  O4 ).  The  temperature  gradient  is 

then 


|^  =  a,=(T,-T,)/(z,-Zj)  (27) 

The  instantaneous  heat  flow  is 

q.  =  -kf  =  -k(T,-Tj/(z,-zJ  (28) 

where  k  is  the  thermal  conductivity  of  the  element  material. 
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The  incremental  increase  in  thermal  energy  at  the  nodes  can  be  obtained  by  integrating  the  heat 
flow  with  respect  to  area  and  time. 

AQi  =  q^AAt + AQj  /  2  (29) 

where  A  s  1 .0  is  the  area,  At  is  the  integration  time  increment,  and  AQe  is  the  internal  energy 
generated  in  the  element  during  the  previous  time  increment. 

After  Equations  26-29  have  been  applied  to  all  elements,  the  updated  temperatures  of  the  nodes 
have  the  form 

=T;+lAQi/MiCpi  (30) 

where  and  T/  are  the  temperature  of  the  nodes  at  times  t+At  and  t,  Z  AQ;  is  the  sum  of  the 

incremental  heat  contributed  by  all  elements  that  contain  node  i,  Mj  is  the  total  mass  of  node  i, 

and  Cpi  is  the  specific  heat  of  node  i. 

The  internal  energy  in  an  element  can  be  used  to  compute  element  pressures.  To  account  for  the 
flow  of  internal  energy  through  the  grid,  the  element  temperature  is  assumed  to  be  the  average  of 
the  nodal  temperatures. 

T  =  (Ti+Tj)/2  (31) 

The  internal  energy  (per  initial  volume)  is  given  by 

E,=(T-T.)p.c,  (32) 

where  To  is  the  initial  temperature  (where  Es  =  0),  po  is  the  initial  density  and  Cp  is  the  specific 
heat  of  the  element  material. 

The  integration  time  increment  must  also  be  bounded  to  ensure  that  the  computations  remain 
stable  for  heat  conduction  (References  6  and  7).  The  heat  conduction  portion  requires 

At  <  pCph^  /  4k  (33) 
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where  h  =  Zj  —  zj  is  the  length  of  the  element  and  the  other  terms  recently  have  been  defined.  This 
is  analogous  to  the  time  increment  restriction  for  wave  propagation  in  Equation  4  of  subsection 
3.1.1.  Unless  h  is  very  small,  the  wave  propagation  restriction  is  much  more  severe  than  the  heat 
conduction  restriction. 

3.1.4  Nonreflective  Boundary  Elements 

Noiueflective  boimdary  elements  have  been  available  for  many  years  (References  8  and  9),  and 
have  been  used  to  provide  nonreflective  damping  for  the  modeling  of  infinite  bodies  such  as  soil 
or  water.  They  are  intended  primarily  to  absorb  elastic  waves,  but  they  are  also  effective  for 
other  applications  (Reference  10).  The  nonreflective  boundary  elements  are  incorporated  as 
single-node,  infmitely-thin,  massless  elements.  The  damping  force  at  node  i  resists  the  velocity 
and  is  expressed  as 

F'=-a„AZi=-p„C3AZi  (1) 

where  the  damping  stress  on  the  boundary  is  =  p„c^,  A  =  1.0  is  the  area  and  Zj  is  the  nodal 
velocity.  The  initial  density  is  p,,  and  the  longitudinal  sound  velocity  is 

c.  =  VK+4G/3)/p.  (2) 

where  Ki  is  the  elastic  bulk  modulus  and  G  is  the  elastic  shear  modulus. 

3.2  ID  CYLINDRICAL  GEOMETRY 

The  element  for  ID  cylindrical  geometry  is  shown  in  Figure  3.  Again,  the  z  coordinate  is  the 
coordinate  of  interest  and  the  element  is  defined  by  nodes  i  and  j,  with  Zi  >  zj.  Only  a  A0  segment 
of  the  element  is  shown  in  Figure  3,  such  that  the  effect  of  the  hoop  stress,  ae,  can  be  illustrated. 

The  nodal  masses  for  the  ID  cylindrical  geometry  are  different  from  the  ID  Cartesian  geometry. 


Mi=aiPoVo 

(1) 

Mj  =  a^PoVo 

(2) 
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Figure  3.  Description  of  the  ID  Cylindricai  Element 


where  po  is  the  initial  density  of  the  material  and  the  initial  volume  is 


V„  =  7i(zf-zj)Ax 


where  Ax  s  1.0  is  a  unit  thickness  normal  to  the  plane  of  the  element  in  Figure  3. 

The  distribution  of  the  total  element  mass  is  determined  by  ai  and  Oj.  The  larger  radius  of  node  i 
requires  a  larger  mass  on  node  i  to  conserve  the  CG  position  of  the  element. 

ttj  =  (2Zi  /3  +  Zj  /3)/(zj  +Zj) 
aj  =  1.0-ai 
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The  strain  rate  in  the  z  direction  is  identical  to  that  of  the  ID  Cartesian  element  in  Equation  4  of 
subsection  3.1.3.  Also,  the  through-thickness  strain  rate  remains  ^  •  The  hoop  strain  rate 

for  this  geometry  is  dependent  on  the  z  velocities. 

ee=z/z  =  (zi  +  Zj)/(zi-i-Zj)  (6) 

This  hoop  strain  rate,  £9 ,  can  be  substituted  for  in  the  ID  Cartesian  strain  rate  equations  in 

subsection  3.1.3.  Also,  the  hoop  stress,  ae,  can  be  substituted  for  ay,  for  the  ID  Cartesian 
stresses. 

The  nodal  forces  depend  on  both  Gz  and  ae- 

F‘  =  -27ia,z-  7cae(zi  -  zj  (7) 

Fj  =  27ta,z-7tae(zj-z.)  (8) 

The  forces  on  nodes  i  and  j  are  not  equal  and  opposite  as  they  are  for  the  ID  Cartesian  geometry, 
and  this  is  due  to  the  hoop  stress,  ae,  that  acts  in  the  same  direction  on  both  nodes.  The  area  for 
the  Gz  stresses  is  taken  at  the  center  of  the  element  at  z  =  ^Zj  +  Zj )  /  2 ,  and  this  provides  equal 

and  opposite  forces  from  the  Gz  stress. 

The  heat  conduction  formulation  and  nonreflective  boundary  elements  are  very  similar  to  those 
of  the  ID  Cartesian  geometry.  The  primary  difference  is  that  the  area.  A,  is  a  function  of  z  for 
the  ID  cylindrical  geometry. 

3.3  ID  SPHERICAL  GEOMETRY 

The  element  for  ID  spherical  geometry  is  similar  to  that  of  ID  cylindrical  geometry  shown  in 
Figure  3.  The  difference  is  that  the  spherical  geometry  does  not  have  a  unit  thickness  normal  to 
the  figure,  but  rather  has  a  A0  and  a  ae  acting  in  the  two  directions  normal  to  the  z  axis. 

The  nodal  masses  are 

(1) 
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Mj=ajp„V„ 


(2) 


where  po  is  the  initial  density  of  the  material  and  the  initial  volume  is 
Vo=47c(zf-z])/3 


(3) 


The  distribution  of  mass  to  the  nodes  (to  conserve  the  CG  position  of  the  element)  is  provided  by 
ai  and  aj. 


Zi  -  T  ZjZf  + 


a.  = 


UJ  ^ 


-(zf-ZjZf-ZiZ]+Z;) 


Oj  =  1.0 -ttj 


(4) 

(5) 


The  nodal  forces  are 

=  -47C<j^z^  -  4jiae  z(zi  -  Zj ) 

FJ  =  47ia^  z^  -  47ca0  z(zi  -  z  ^ 

Again,  the  area  for  the  Oz  stresses  is  taken  at  z  =  (z^  +  Zj)  /  2 .  Other  features  of  this  geometry  are 
analogous  to  the  ID  Cartesian  and  cylindrical  geometries. 


3.4  2D  AXISYMMETRIC  GEOMETRY 

A  description  of  2D  axisymmetric  geometry  is  shown  in  Figure  4.  This  figure  shows  triangular 
elements,  but  quad  elements,  bending  shell  elements,  membrane  shell  elements,  and 
nonreflective  boundary  elements  are  also  available.  The  early  developments  of  this  work  are 
reported  in  References  1, 3,  and  1 1. 
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3.4.1  Equations  of  Motion 


The  radial  (x)  and  axial  (z)  accelerations,  velocities,  and  positions  of  node  i  are  determined  as 
follows: 


x‘  =  I;VM< 

(1) 

z‘=F'/Mi 

(2) 

x'*  =x|"  +  xjAt 

(3) 

z|^=z‘-+z|At 

(4) 

=  X*  +  x*"^  At 

(5) 

z‘^^‘=z|  +  zrAt 

(6) 

Equations  1  and  2  provide  the  radial  and  axial  accelerations  at  time  =  t.  and  F^'  are  the  net  x 
and  z  direction  forces  on  node  i.  Looking  at  the  lower  portion  of  Figure  4,  these  forces  would 
come  from  the  four  elements  that  contain  node  i.  The  forces  are  obtained  from  the  element 
computations  in  the  previous  cycle  of  integration.  Mj  is  the  total  mass  at  node  i,  and  it  contains 

a  fraction  of  the  mass  of  all  elements  that  contain  node  i.  Although  Figure  4  shows  only 
triangular  elements,  a  node  can  have  various  element  types  attached  to  it.  This  means  it  is 
possible  for  an  individual  node  to  have  masses  and  forces  from  different  element  types.  The 
specific  algorithms  for  masses  and  forces  are  provided  in  subsections  3.4.4-3.4.8. 

The  updated  velocities  in  Equations  3  and  4  are  constant  for  the  interval  between  time  =  t  a^ 
time  =  t+At.  The  constant  velocities  for  the  previous  time  increment  are  xj"  and  z*" ,  and  At  is 

the  average  of  the  two  integration  time  increments  about  time  =  t. 

The  integration  time  increment  is  limited  to 


At  =  C.[h/(V?+Vg'+cO. 


(7) 
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where  =  CqQ/p,  h  is  a  characteristic  length  of  the  element  (described  later  for  each  element 
type),  and  Cs  is  the  sound  velocity  of  the  material  (References  1  and  5).  Q  is  the  artificial 
viscosity  and  Cq  is  a  constant  for  the  artificial  viscosity.  These  are  described  in  Section  6.  The 
Courant  sound  speed  firaction,  Q,  must  be  less  than  unity  (Ct  <  1.0)  to  ensure  numerical  stability. 
Ct  =  0.9  is  a  typical  value. 

For  2D  bending  shell  elements  there  is  also  a  rotational  degree  of  fi'eedom  on  each  node.  The 
rotational  velocity  of  node  i  is  updated  in  a  similar  manner  as  the  translational  velocities  in 
Equations  1-4.  The  updated  acceleration  and  velocity  are  as  follows: 

e'=M;/iy  (8) 

e'+=e‘-+e'At  (9) 

where  MJ,  is  the  net  moment  on  node  i,  ly  is  the  rotational  mass  moment  of  inertia,  and  At  is  the 
average  integration  time  increment.  The  rotational  displacement,  0i,  is  not  required. 

3.4.2  Sliding  Interfaces 

The  2D  sliding  algorithm  allows  the  two  surfaces  (master  and  slave)  to  be  determined 
automatically,  or  to  be  input  by  the  user.  This  subsection  describes  the  automated  approach  in 
detail.  It  is  similar  to  that  reported  in  Reference  12,  but  it  has  been  improved  to  eliminate  the 
order  dependence.  There  are  three  primary  steps  in  the  automated  sliding  interface  algorithm;  the 
interface  determination  algorithm,  the  searching  algorithm  and  the  contact  algorithm. 

Interface  Determination  Algorithm  —  After  the  finite  element  grid  has  been  assembled  (using 
three-node  triangles  or  four-node  quads),  it  is  necessary  to  determine  the  surface  nodes  and  the 
surface  line  segments.  The  surface  nodes  are  designated  as  slave  nodes,  and  the  surface  line 
segments  are  designated  as  master  segments  (composed  of  two  master  nodes).  The  left  portion 
of  Figure  5  shows  a  simple  illustration  for  an  impact  problem.  The  assemblage  of  master 
segments  forms  a  membrane  around  the  outer  surfaces  of  the  bodies  through  which  the  slave 
nodes  cannot  pass.  As  the  slave  nodes  contact  the  master  segments,  momentum  is  transferred 
from  the  slave  nodes  to  the  master  nodes.  Note  that  every  surface  node  is  both  a  slave  node  and  a 
master  segment  node.  This  is  different  from  the  earlier  sliding  algorithms,  and  therefore  requires 
special  treatment. 
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The  right  portion  of  Figure  5  shows  how  the  number  of  slave  nodes  and  master  segments  can  be 
significantly  reduced  by  the  user,  if  the  user  has  some  knowledge  about  the  region  of  interaction. 
Although  this  is  not  required,  it  can  often  be  a  simple  way  to  reduce  the  number  of  slave  nodes 
and  master  segments  in  a  problem.  This,  in  turn,  leads  to  reduced  CPU  time. 


'  All  surface  segment  lines  (2  nodes) 
designated  as  master  segments 

T19616 


Figure  5.  Determination  of  the  2D  External  Interfaces 

The  master  segments  can  be  obtained  by  comparing  every  line  segment  of  every  element  with 
every  line  segment  of  every  other  element.  If  there  is  no  matching  line  segment  (with  the  same 
two  nodes),  then  that  segment  is  on  the  surface  and  is  designed  a  master  segment.  Likewise, 
every  node  that  forms  a  master  segment  is  also  a  slave  node.  The  final  description  of  the 
interfaces  consists  of  a  list  of  slave  nodes  and  a  list  of  master  segments  that  are  defined  by  two 
master  nodes.  The  pairs  of  master  nodes  (Mi,  M2)  are  ordered  such  that  the  exterior  surface  is  to 
the  left,  and  the  interior  is  to  the  right,  when  viewing  node  Mi  from  node  M2. 


Searching  Algorithm  —  After  the  nodal  equations  of  motion  (velocities  and  positions)  have 
been  updated,  it  is  necessary  to  check  all  of  the  slave  nodes  to  determine  if  they  have  interacted 
with  any  of  the  master  segments  that  form  the  master  surface.  A  summary  of  the  searching 
algorithm  is  shown  in  Figure  6.  The  master  nodes  are  Mi  and  M2,  and  the  slave  node  is  Ns.  The 
first  check  for  each  slave  node  is  to  determine  which  master  segment  (or  segments)  are 
candidates  for  interaction.  A  bucket  sort  algorithm  (Reference  13)  is  used  to  limit  the  number  of 
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master  segments  that  must  be  considered  for  each  slave  node.  If  one  of  the  two  master  nodes  on 
a  segment  is  also  the  slave  node,  then  that  master  segment  is  no  longer  considered  for  that  slave 
node.  The  box  test  in  Figure  6  indicates  that  slave  node  Ns  can  be  associated  with  master 
segment  M1-M2  only  if  it  is  contained  within  a  rectangular  box  that  extends  a  distance,  5ref, 
beyond  the  master  segment,  M1-M2. 


Figure  6.  Description  of  the  2D  Sliding  Interface  Searching  Algorithm 


The  reference  distance  is  simply 


(1) 


where  Vref  is  the  maximum  relative  velocity  (between  a  slave  node  and  a  master  segment)  that 
can  be  experienced  during  the  cotirse  of  the  problem,  and  At  is  the  integration  time  increment. 
Vref  can  be  specified  by  the  user  or  it  can  be  computed  by  the  preprocessor  based  on  the 
maxirmim  relative  velocity  for  the  initial  condition.  It  can  generally  be  set  to  about  1 .5  times  the 
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maximum  relative  velocity  in  the  initial  condition.  For  explosives  it  can  generally  be  set  to  the 
detonation  velocity. 

If  the  slave  node  passes  the  box  test,  it  is  then  subjected  to  the  cross-over  test.  This  test  is 
performed  to  determine  if  the  slave  node  has  crossed  over  the  master  segment.  The  distance 
between  the  slave  node  and  the  master  segment  is  given  by 

5  =  -(AXj  +  BZj  +  c)  (2) 

where  A  =  (z2-Zj)/f ,  B  =  (x, -Xj)/^,  C  =  (x2Z, -XiZj)/^ ,  and  ^  =  -y/(x,  -X2)  -i-(z, -Z2)  . 

The  X  and  z  coordinates  of  master  nodes  Mi  and  M2  are  xi,  X2,  zi,  Z2,  and  the  coordinates  of  slave 
node  Ns  are  Xs  and  Zs.  A  and  B  are  the  direction  cosines  normal  to  the  master  segment  and  i  is 
the  length  of  the  segment.  If  0  <  5  <  5ref,  then  the  slave  node  has  touched  or  crossed  over  the 
master  segment  line  (or  the  extension  of  the  line),  and  it  remains  a  candidate  for  interaction  with 
that  master  segment.  If  the  slave  node  fails  either  the  box  test  or  the  cross-over  test,  then  it 
cannot  interact  with  that  master  segment. 

If  the  slave  node  passes  both  the  box  test  and  the  cross-over  test,  then  it  must  be  determined  if  it 
falls  within  the  normal  projection  of  master  segment  Mi— M2  (Region  A),  or  if  it  falls  within  the 
extended  region  of  the  master  segment  (Region  B).  Even  after  a  slave  node  has  found  a  master 
segment  (by  passing  both  the  box  test  and  the  cross-over  test),  the  search  must  continue  for  the 
remaining  master  segments,  because  it  is  possible  for  a  slave  node  to  find  multiple  candidate 
master  segments. 

It  is  also  necessary  to  check  if  the  slave  node  has  not  crossed  over  the  master  segment,  but  is 
within  Region  C  (-5ref  <  5  <  0),  as  shown  in  Figure  6.  This  condition  can  be  used  to  offset  a 
Region  A  condition  for  an  adjacent  segment,  as  will  be  shown  later. 

After  the  slave  node  has  searched  for  all  the  master  segments,  there  are  the  following 
possibilities: 

•  If  the  slave  node  is  within  Region  A  for  only  one  master  segment  and  is  not  within 
Region  C  for  an  adjacent  segment,  then  it  interacts  with  master  segment  Mb-Mc,  as 
shown  in  the  lower  portion  of  Figure  6. 
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•  If  a  slave  node  is  within  Region  A  for  only  one  master  segment  and  is  within  Region  C 
for  an  adjacent  segment,  then  it  may  or  may  not  be  in  contact  with  the  master  surface. 
This  condition  occurs  in  the  shaded  area  near  node  Md,  as  shown  in  the  lower  portion  of 
Figure  6.  If  the  slave  node  passed  through  segment  Mc-Md,  then  it  is  designated  to 
interact  with  that  segment.  If  it  did  not  pass  through  segment  Mc-Md,  then  it  is  not  in 
contact  with  the  master  surface. 

•  If  the  slave  node  is  within  Region  A  for  two  adjacent  master  segments,  then  it  interacts 
with  whichever  segment  it  crossed  to  arrive  at  its  current  position.  This  possibility  is 
illustrated  in  the  shaded  area  near  node  Me,  in  the  lower  portion  of  Figure  6.  If  it  passed 
through  line  Mb-Mc  to  reach  its  position  in  the  shaded  area,  then  it  interacts  with  master 
segment  Mb-Mc.  Likewise,  if  it  passed  through  line  Mc-Md,  then  it  interacts  with  that 
segment. 

•  If  the  slave  node  is  within  Region  B  for  only  one  segment,  then  it  does  not  interact  with 
any  master  segment.  This  condition  can  exist  at  a  convex  comer  on  the  master  surface. 
The  node  can  be  crossed  over  the  extension  of  a  master  segment  (Region  B),  but  it  may 
not  be  crossed  over  the  adjacent  master  segment. 

•  If  the  slave  node  is  within  Region  B  for  two  adjacent  master  segments,  then  it  interacts 
only  with  the  single  node  that  is  common  for  those  two  segments.  This  possibility  is 
illustrated  in  the  shaded  area  near  node  Mb,  in  the  lower  portion  of  Figure  6. 

•  Another  situation  can  exist  that  is  not  illustrated  in  Figure  6.  Under  some 
circumstances,  where  three  or  more  distinct  bodies  meet  at  a  comer,  it  is  possible  for  a 
slave  node  to  pass  through  two  or  more  master  segments  that  are  not  attached  to  one 
another.  When  this  occurs,  the  slave  node  is  associated  with  the  master  segment  that 
has  the  greatest  crossover.  If  the  crossover  with  other  master  segments  is  significant 
(6  >  0.1  6ref)  then  the  complete  searching  and  contact  procedure  is  repeated  for  that 
cycle,  with  the  result  being  that  the  secondary  crossover  now  becomes  the  primary 
crossover  and  is  adjusted  accordingly. 
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After  the  searching  is  completed,  there  are  three  possibilities  for  each  slave  node. 


«  The  slave  node  has  found  a  master  segment  M1-M2  with  which  it  must  interact. 


.  The  slave  node  has  foimd  a  single  master  node  with  which  it  must  interact. 

•  The  slave  node  has  not  crossed  over  any  portion  of  the  master  surface,  and  therefore  has 
no  interaction  with  the  master  surface. 

The  searching  algorithm  is  applied  every  cycle.  In  some  instances,  the  CPU  searching  time  can 
be  reduced  by  performing  the  bucket  sort  algorithm  every  N  cycles  instead  of  every  cycle.  This 
requires  the  buckets  to  be  expanded  in  size  such  that  all  the  candidate  master  segments  are 
available  to  the  slave  node  for  all  N  cycles.  This  is  accomplished  by  using  the  reference  velocity 
concept  of  Equation  1.  This  approach  provides  greater  CPU  savings  for  lower  reference 
velocities  because  the  buckets  do  not  need  to  be  expanded  as  much  as  they  are  for  higher 
reference  velocities.  This  approach  cannot  be  used  with  erosion,  however,  because  the  lists  of 
slave  nodes  and  master  segments  ehange  from  cycle  to  cycle  during  erosion. 

Contact  Algorithm  —  The  contact  algorithm  adjusts  the  positions  and  velocities  of  the  slave 
nodes  and  the  master  surface  nodes.  This  is  done  after  the  equations  of  motion  have  been 
updated  in  the  standard  manner,  and  after  the  searching  has  been  completed.  Conservation  of 
momentum  can  be  attained  by  considering  each  slave  node  to  interact  with  its  master  segment  (or 
single  master  node)  in  a  single  pass  (without  subsequent  interactions).  However,  subsequent 
interactions  generally  will  be  required  to  have  each  slave  node  placed  exactly  on  the  master 
surface,  and  to  have  the  normal  velocity  of  the  slave  node  equal  to  the  normal  velocity  of  the 
master  sxirface  at  the  slave  node  position. 

A  summary  of  the  contact  algorithm  is  shown  in  Figure  7.  The  general  approach  is  to  have  each 
slave  node  interact  with  its  master  segment  (Mi  -  M2)  or  master  node  (if  it  is  in  a  concave 
comer).  The  velocity  and  position  changes  for  the  master  and  slave  nodes  are  performed  during 
the  course  of  several  iterations.  Multiple  iterations  are  required  to  achieve  a  good  velocity  and 
position  match  of  the  slave  node  with  the  master  segment.  For  a  single  slave  node  on  a  single 
master  segment,  it  is  possible  to  obtain  an  exact  match  in  one  iteration.  However,  when 
subsequent  slave  nodes  interact  with  the  same  master  segment,  or  an  adjacent  segment,  then  the 
match  of  the  previous  slave  node  is  altered. 
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Intermediate 


Three  adjusted  normal  velocities 
(Vi,  Vj,  V3)  require  three  equations: 

•  Conserve  normal  momentum 

•  Conserve  angular  momentum 

•  Equate  normal  slave  node  velocity  to  master 
segment  velocity  at  slave  node  position 


Figure  7.  Description  of  the  2D  Sliding  Interface  Contact  Algorithm 


Furthermore,  only  a  portion  of  the  required  velocity  and  position  changes  are  made  during  each 
iteration.  If  the  entire  position  movement  would  be  made  during  the  first  iteration,  then  it  would 
be  possible  for  some  slave  nodes  to  be  separated  fi*om  the  master  segment  after  the  first  iteration. 
Because  the  contact  algorithm  is  exercised  only  when  there  is  interference,  if  the  slave  node  is 
free  from  the  surface,  then  it  will  not  be  considered  for  subsequent  iterations. 

An  alternative  approach  would  be  to  consider  the  slave  node  (and  the  corresponding  master 
segment)  to  be  adjusted  for  all  subsequent  iterations,  if  it  was  in  contact  for  the  first  iteration. 
The  problem  with  this  approach  is  that  this  would  not  allow  a  slave  node  to  be  released  from  the 
master  surface  if  it  is  in  contact  during  the  first  iteration.  For  some  complex  contact  interfaces,  a 
slave  node  needs  to  be  free  to  separate  from  the  master  surface,  even  if  it  was  in  contact  during 
the  first  iteration. 


The  following  algorithm  considers  a  contact/sliding  interface  with  no  friction,  although  frictional 
effects  can  be  added  straightforwardly  using  the  general  approach  described  in  Reference  14,  and 
presented  in  subsection  3.8.2.  The  basic  approach  is  to  consider  each  slave  node  interacting  with 
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its  master  segment.  The  normal  velocities  of  the  slave  node  and  the  two  master  nodes  are  V," , 

V," ,  and  Vj" ,  where  the  superscript  n  indicates  the  iteration  number.  The  objective  of  the 
algorithm  is  to  adjust  these  normal  velocities  by  conserving  linear  momentum,  conserving 
angular  momentum,  and  equating  the  normal  velocity  of  the  slave  node  to  the  normal  velocity  of 
the  master  segment  at  the  slave  node  location. 

The  relationship  between  the  normal  velocity  changes  (based  on  conservation  of  linear  and 
angular  momentum)  is  as  follows: 

A\?=-R;M.AV."/M„j  W 

where  AV"  and  AVj"  are  the  normal  velocity  changes  to  master  nodes  Mi  and  M2  during  the  n‘*’ 

iteration;  Ri  and  R2  are  the  fractions  of  momentum  transferred  from  the  slave  node  to  the  master 
nodes;  Ms,  Mmi,  and  Mm2  are  the  masses  at  the  slave  node  and  the  two  master  nodes;  and  AV^  is 

th 

the  normal  velocity  change  to  Ae  slave  node  during  the  n  iteration. 


Ri  and  R2  are  simply: 

R,  =  A/(^sn.-X2)'+(Zsm-Z2  )'  '  ^ 

R2=1.0-Ri 


(5) 

(6) 


where  Xsm  =  Xs  +  A5  and  Zsm  =  Zs  +  B5  are  the  coordinates  of  the  slave  node  when  projected  back 
to  the  master  segment. 


Now,  if  the  normal  velocity  of  the  slave  node  is  equated  to  the  normal  velocity  of  the  master 
segment  (at  the  slave  node  location),  and  the  relationships  of  Equations  3  and  4  are  substituted, 
the  normal  velocity  change  to  the  slave  node  is: 


Avr*  =• 


cx(v:-v.‘) 


(i  +  R^M^/Mmi+RjM, 


/M„j) 


(7) 
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where  a  =  1 .0  for  an  exact  match  and  V“  =  Y"  +  (^i"  -  V2’' )  is  the  master  normal  velocity  at 

the  projected  slave  node  location. 

In  Equation  7,  it  is  possible  to  substitute 

8/At=(v;-v.-)  (8) 


which  gives 


Avr‘ 


_ a(S  /  At) _ 

(l  +  R;M./M„,+RjM./M„,) 


(9) 


where  5  is  the  normal  deflection  from  Equation  2  and  At  is  the  integration  time  increment.  The 
use  of  Equations  8  and  9,  instead  of  Equation  7,  allows  the  algorithm  to  be  self-correcting  when 
the  previous  cycle  does  not  provide  a  perfect  velocity  and  position  match. 

As  was  noted  previously,  it  is  necessary  to  perform  only  a  fraction  of  the  velocity  and  position 
change  during  each  iteration.  This  generally  keeps  the  slave  node  in  a  crossed-over  position  for 
intermediate  iterations,  as  shown  in  Figure  7,  such  that  subsequent  iterations  can  gradually  move 
it  toward  the  master  segment.  However,  if  the  slave  node  moves  free  during  an  intermediate 
iteration,  then  it  is  not  adjusted  during  the  next  iteration.  The  expression  used  for  a  is  taken  as 

a  =  aia2  (10) 


The  first  factor  is  based  on  the  iteration  nmnber. 


a,  = 


1 


'  VN-n+1 


(11) 


where  N  is  the  total  number  of  iterations  and  n  is  the  current  iteration  number.  Note  that  ai  =  1 
when  n  =  N. 


The  second  factor  is  based  on  the  number  of  slave  nodes  that  act  on  node  Ns,  when  it  (Ns)  is  a 
master  node.  This  is  expressed  as 
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IR„+IR. 


(12) 


where  ZRs  is  the  sum  of  Ri  and  R2  (from  Equations  5  and  6)  that  node  s  receives  from  other  slave 
nodes  when  it  (Ns)  is  acting  as  a  master  node.  ZRm  is  a  similar  sum  on  the  master  surface.  This 
is  required  to  keep  nodes  from  being  adjusted  too  much  when  they  are  both  a  slave  and  a  master. 

Note  that  if  slave  node  Ns  is  never  a  master  node,  then  ERs  =  0  and  aa  =  1 .0.  This  means  that  the 
basic  algorithm  can  be  used  for  the  following  three  categories: 

•  Automatic 

•  Symmetric 

•  Nonsymmetric. 

For  automatic  sliding  all  of  the  surfaces  are  automatically  designated  as  slave  nodes  and  master 
surfaces  and  the  user  does  not  need  to  specify  the  interfaces.  This  is  the  recommended  option. 

For  symmetric  sliding  both  sides  of  the  interface  are  designated  as  master  surfaces  and  slave 
nodes.  This  gives  a  symmetric  interface  that  is  not  order  dependent  and  is  not  dependent  on 
which  side  is  master  or  slave.  This  option  provides  the  same  results  as  the  automatic  option  if  all 
of  the  interacting  surfaces  are  designated  as  master  surfaces  and  slave  nodes. 

The  nonsymmetric  option  is  the  traditional  option  where  one  side  is  the  master  surface  and  the 
other  side  is  slave  nodes.  Generally,  the  master  surface  should  have  the  denser  and  stronger 
material,  an  equal  or  greater  node  spacing  than  the  slave  nodes,  and  should  not  have  a  convex 
surface  toward  the  slave  nodes.  This  is  the  condition  that  exists  for  an  SPH  node,  or  a  standard 
node  whose  elements  have  eroded,  that  interacts  with  a  master  surface. 

After  the  slave  velocity  change  has  been  determined  from  Equation  9,  then  the  x  and  z  velocity 
changes  to  the  slave  node  are  then  transferred  back  into  the  system  coordinates 


Au^^’  =AAV,"^' 

(13) 

Av^’  =  BAV”^' 

(14) 
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The  velocity  changes  Av"'"*,  Au"'"\  Av"""')  to  the  master  nodes  are  performed  in  a 

manner  similar  to  that  of  Equations  13  and  14,  where  AV" ,  and  AV2"  are  obtained  by  substituting 
AVj”  into  Equations  3  and  4. 

After  the  velocity  changes  have  been  made  to  all  nodes  on  the  sliding  interface,  then  the  positions 
of  the  nodes  must  be  updated  in  a  manner  that  is  consistent  with  the  velocity  changes.  The 
updated  nodal  positions  for  all  nodes  on  the  sliding  interface  are 

x."^‘=x"+lAu;‘^‘At  (15) 

=  z"  +  S  Av|''^' At  (16) 

where  S  Au"^'  and  S  Avf^’  are  the  sum  of  all  master  and  slave  velocity  changes  made  to  node  i. 
During  this  process  there  is  no  arbitrary  movement  of  nodes,  but  rather  total  consistency  between 
the  velocity  changes  and  position  changes.  Because  the  velocity  and  position  changes  are  made 
after  all  the  interface  nodes  are  processed,  there  is  no  order  dependence. 

The  application  of  the  preceding  algorithm,  for  all  slave  nodes,  completes  a  single  iteration. 
Subsequent  iterations  are  performed  in  a  similar  manner.  The  number  of  iterations  has  little 
effect  on  the  CPU  time  because  the  iterations  require  much  less  time  than  the  searching. 

3.4.3  Erosion 

For  many  problems  it  is  necessary  to  allow  highly  distorted  material  to  be  eroded  away.  This  is 
done  to  simulate  the  actual  phenomenon  that  occurs  for  erosion  problems,  and  to  delete  the 
highly  distorted  elements  fi-om  the  computation.  This  algorithm  is  reported  in  References  15  and 
16.  It  is  only  available  for  triangular  elements. 

The  algorithm  which  follows  should  be  applied  only  to  problems  where  erosion  is  the  primary 
mode  of  penetration.  Figure  8  shows  a  master  surface  defined  as  a  set  of  consecutively  linked 
nodes  (Ml,  M2  ...  M15,  M16).  The  slave  nodes  are  always  to  the  left  of  the  master  surface  when 
moving  firom  the  beginning  of  the  master  surface  (node  Ml)  to  the  end  of  the  master  surface 
(node  Ml 6).  The  slave  nodes  can  be  randomly  ordered  and  are  not  allowed  to  penetrate  through 
the  master  sxuface. 
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The  master  surface  is  allowed  to  erode  by  totally  failing  triangular  elements  which  have  one,  two 
or  three  sides  on  the  master  surface.  Totally  failed  elements  cannot  develop  any  stresses  or 
pressures;  they  essentially  disappear  except  that  mass  is  retained  at  the  nodes. 

The  procedure  begins  by  putting  a  flag  on  the  nodes  which  are  on  the  master  surface.  Then  the 
three  nodes  of  each  element  are  examined  to  determine  if  the  element  is  on  the  master  surface.  If 
two  or  three  of  the  nodes  on  the  specific  element  are  flagged,  then  that  element  is  on  the  master 
surface  and  it  may  be  totally  failed  if  it  meets  one  of  the  failure  criteria. 

The  five  elements  in  Figure  8  (elements  A,  B,  C,  D,  E)  have  two  or  three  consecutive  nodes  on 
the  master  surface,  and  have  one,  two  or  three  sides  on  the  master  surface.  The  following 
describes  the  criteria  necessary  to  totally  fail  these  elements  and  to  redefine  the  master  surface 
and  the  slave  nodes. 

.  Element  A  has  one  side  (M1-M2)  on  the  master  surface  and  one  side  (Ml-Nl)  on  the 
rotational  axis  of  symmetry.  If  the  equivalent  plastic  strain  or  volumetric  strain  of 
element  A  exceeds  a  specified  eroding  strain,  the  element  is  totally  failed.  The  previous 
master  surface  (Ml,  M2,  M3  ...)  is  then  replaced  by  an  updated  master  surface  (Nl, 

M2,  M3  ...).  Also,  node  Ml  is  designated  as  a  slave  node  to  the  updated  master 
surface. 

.  Element  B  has  one  side  (M3-M4)  on  the  master  surface.  If  the  equivalent  plastic  strain 
or  volumetric  strain  of  element  B  exceeds  a  specified  eroding  strain,  the  element  is 
totally  failed.  The  previous  master  surface  (. . .  M2,  M3,  M4,  M5  . . .)  is  then  replaced 
by  an  updated  master  surface  (. . .  M2,  M3,  N2,  M4,  M5  . . .)  There  are  no  new  slave 
nodes  to  be  designated  for  this  case. 

.  Element  C  has  two  sides  (M5-M6,  M6-M7)  on  the  master  surface.  If  the  equivalent 
plastic  strain  or  volumetric  strain  of  element  C  exceeds  a  specified  eroding  strain,  the 
element  is  totally  failed.  The  previous  master  surface  (. . .  M4,  M5,  M6,  M7,  M8  . . .)  is 
then  replaced  by  an  updated  master  surface  (. . .  M4,  M5,  M7,  M8  . . .).  Also,  node  M6 
is  designated  as  a  slave  node  to  the  updated  master  surface. 
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Figure  8.  Description  of  the  2D  Erosion  Algorithm 

.  Element  D  also  has  two  sides  (M8-M9,  M9-M10)  on  the  master  surface.  If  the  angle 
a,  defined  by  the  two  master  sides,  is  less  than  a  specified  angle,  the  element  is  totally 
failed.  The  previous  master  surface  (. . .  M7,  M8,  M9,  MIO,  Ml  1  . . .)  is  then  replaced 
by  an  updated  master  surface  (. . .  M7,  M8,  MIO,  Ml  1  . . .).  Also,  node  M9  is  designated 
as  a  slave  node  to  the  updated  master  surface.  This  criterion  tends  to  give  a  smoothed 
master  surface,  by  eliminating  intruding  elements  which  have  not  exceeded  the  previous 
eroding  strain  criterion. 

.  Element  E  has  three  sides  (Ml  1-M12,  M12-M13,  M13-M14)  on  the  master  surface. 
This  conditions  exists  because  node  Ml  1  is  identical  to  node  M14.  Under  these 
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conditions,  element  E  is  always  totally  failed.  The  previous  master  surface  (. . .  MIO, 
Ml  1,  M12,  M13,  M14,  M15  .. .)  is  replaced  by  an  updated  master  surface  (...  MIO, 
MU,  M15  ...).  Also,  nodes  M12  and  M13  are  designated  as  slave  nodes  to  the  updated 
master  surface.  This  criterion  also  tends  to  give  a  smoothed  master  surface. 

For  axisymmetric  geometry,  when  the  master  nodes  are  deleted  from  the  master  surface  and 
designated  as  slave  nodes,  the  outward  radial  velocity  of  the  slave  node  is  redefined  to  be  a 
fraction  of  the  axial  velocity  of  that  node.  This  is  done  to  move  the  slave  nodes  initially  on  the 
axis  of  symmetry,  away  from  the  axis,  as  would  be  expected  during  an  actual  eroding  penetration 
process. 

To  ensure  that  there  is  no  cross-over  of  material  between  the  two  surfaces  (such  as  a  projectile 
and  a  target),  a  two-step  approach  (two  separate  sliding  interfaces)  is  used  with  the  standard  (not 
automatic)  sliding  algorithm.  For  the  first  interface,  all  of  the  projectile  nodes  and  the  eroded 
target  nodes  are  slave  to  the  master  surface  in  the  target.  For  the  second  interface,  all  of  the 
nodes  in  the  center  of  the  target,  and  the  eroded  projectile  nodes,  are  slave  to  the  master  surface 
on  the  projectile.  The  eroded  nodes  are  therefore  restrained  from  passing  through  either  the 
projectile  or  the  target.  The  automatic  sliding  algorithm  handles  both  surfaces  automatically. 

3.4.4  Triangular  Elements 

The  axisymmetric  triangular  element  algorithm  follows  the  same  sequence  as  the  ID  Cartesian 
element.  The  element  is  geometrically  defined  by  nodes  i,  j,  and  m  (counterclockwise)  as  shown 
in  Figure  4.  This  algorithm  is  reported  in  References  1, 3,  and  11.  The  coordinates  of  node  i  are 
designated  Xi  and  Zi,  and  the  corresponding  velocities  are  designated  Uj  and  Vj  (which  are 
identical  to  Xj  and  Zj  used  previously). 

The  algorithm  that  follows  is  for  a  single  triangular  element.  The  arrangement  of  the  elements, 
however,  can  have  a  significant  effect  on  the  computed  response  (References  16, 17,  and  18). 

The  crossed  triangle  arrangement  (four  triangles  within  a  quad),  in  the  lower  portion  of  Figure  4, 
generally  provides  good  accuracy.  A  slashed  triangle  arrangement  (two  triangles  within  a  quad) 
can  produce  significant  inaccuracies.  There  are  also  some  mixed  (average  volumetric  strain,  or 
average  pressure)  algorithms  that  provide  improved  acciuucy  for  various  arrangements  of 
elements  (Reference  16). 
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The  mass  at  node  i,  for  an  individual  element,  is 


Mi=aiP„V„ 


(1) 


where  po  is  the  initial  density  of  the  material  and  Vo  is  the  initial  volume  of  the  element.  The 
fraction  of  the  element  mass  that  is  assigned  to  node  i  is 


1  irx,.^ 


4  "^12 


(2) 


where  the  average  radial  coordinate  is 


x  =  (xi+Xj  +  x^)/3  (3) 

The  initial  volume  is 

V„  =  27i5cA„  (4) 


where  Ao  is  the  initial  cross-sectional  area. 

When  element  A  is  incorporated  into  an  assemblage  of  elements,  as  shown  in  Figure  4,  then  the 
total  mass  at  node  i  contains  individual  element  masses  from  all  elements  that  contain  that  node. 
The  total  mass  at  node  i  is 

M,=SM|  (5) 


The  volumetric  strain  and  strain  rate  are  obtained  in  the  same  manner  as  used  for  the  ID 
Cartesian  element. 

E,=V/V,-1  (6) 

e,  =  (e""-e;)/At  (7) 

In  Equation  6,  V  and  Vo  are  the  current  and  initial  element  volumes;  and  in  Equation  7,  and 
are  the  volumetric  strains  at  time  =  t  +  At  and  time  =  t.  The  volumetric  strains  and  strain  rates 
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are  based  on  the  initial  configuration,  as  opposed  to  the  shear  and  deviator  strain  rates  which  are 
based  on  the  current  configuration. 

The  strain  rates  are  obtained  fi-om  the  current  geometry  of  the  element  and  the  velocities  of  the 
nodes.  If  it  is  assumed  that  the  lines  connecting  the  nodes  remain  straight,  the  displacements  and 
velocities  within  each  element  must  vary  in  a  linear  manner.  Then  the  velocities  within  the 
element  can  be  expressed  as 

u  =  a,  H-ttjX+ajZ 

y  =  a^+aiX+a^z 


where  ai ...  as  are  geometry  and  velocity-dependent  constants.  It  is  possible  to  solve  for  ai,  az, 
as  by  substituting  the  radial  velocities  and  coordinates  of  nodes  i,  j,  m  into  Equation  8.  This 
gives  three  equations  and  three  unknowns  such  that  Equation  8  can  be  expressed  in  terms  of  the 
element  geometry  and  nodal  velocities. 

u  =  — [(a;  -t-biX-i-Ciz)ui  -i-(aj  +bjX+Cjz)uj  -i-(a„  +b„x+c„z)u^  j  (10) 

2A  ^ 

where  a^  =  XjZ„  -  x^z^ ,  b-  z^  -  z„ ,  c,  =  x„  -  x^ ,  and  A  is  the  cross-sectional  area  of  the 

element  in  the  x-z  plane.  The  axial  velocities  are  identical  to  Equation  10  except  the  radial 
velocities  are  replaced  by  the  axial  velocities.  Another  geometry  constant  for  triangular  elements 
is  the  miTiirmitn  altitude  of  the  triangle,  h.  This  is  used  for  the  integration  time  increment  in 
Equation  7  of  subsection  3.4.1. 

After  the  velocities  are  obtained,  it  is  possible  to  determine  the  normal  strain  rates  » 

the  shear  strain  rate  (Yxz)  Oie  localized  rotational  spin  rate  of  the  element  in  the  x-z  plane 

(®xz)- 


1  r  .  .  .1 


dz  2A 


(11) 


(12) 
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(13) 


ee=|  =  (ui+Uj+u,)/(xi+Xj+x,) 

3u  9v 


(14) 


_  1  ^  av 
2  tax  az; 


(15) 


It  can  be  seen  that  Equations  1 1, 12, 14,  and  15  are  derivatives  of  linear  functions  and  are 
therefore  constant  within  the  element.  Equation  13  involves  averages  of  the  nodal  velocities  and 
the  three  radii  (h  and  x) ,  so  it  is  not  necessarily  constant. 

The  other  two  shear  strain  rate  components  are  zero  =  Yze  =  o)  for  axisymmetric  geometry 

with  no  spin  about  the  z  axis  of  rotation.  These  become  nonzero  when  spin  is  included,  and  this 
is  described  in  subsection  3.5. 

The  equivalent  strain  rate  (e)  and  deviator  strain  rates  (e,, 6^,69)  are  determined  in  an  identical 
manner  as  shown  previously  for  ID  geometry. 


^ave 

where  £  +  Ee )  /  3 .  Note  that  the  sum  of  the  deviator  strain  rates  is  ex+e2+eQ=0. 

The  stresses  in  the  elements  are  determined  from  the  strains,  strain  rates,  temperatures,  pressures, 
internal  energies,  and  material  constants.  This  is  identical  to  the  ID  algorithm,  except  that  a 
shear  stress  is  also  present.  The  three  normal  stresses  are  generally  expressed  as 


A  t19616.doc 


41 


<5.=s.-(P  +  Q) 

(20) 

5.=S,-(P  +  Q) 

(21) 

<5e  =  s,-(P  +  Q) 

(22) 

where  Sx,  Sz,  and  se  are  the  normal  deviator  stresses,  P  is  the  hydrostatic  pressure  and  Q  is  the 
artificial  viscosity.  These  are  described  in  detail  in  Section  6.  The  nonzero  shear  stress  is  Xxz  and 
the  two  zero  shear  stresses  for  this  geometry  are  Xxe  and  Xze. 

Trial  values  of  the  deviator  stresses  at  time  =  t  +  At  are 

s‘^^'=s‘+2Ge,At-2x>^At 

(23) 

s^^'  =  s^  +  2Ge,  At + 2x>^At 

(24) 

Sg'^^  =  Sg  +  2Geg  At 

(25) 

In  Equation  23  the  first  term  (s^)  is  the  radial  stress  at  the  previous  time  and  the  second  term 
(2Ge^At)  is  the  incremental  stress  due  to  the  incremental  strain  (e^At)  during  that  time 
increment,  where  G  is  the  elastic  shear  modulus.  The  third  term  (2T^®,„.At)  is  due  to  shear 
stresses  from  the  previous  time  increment,  which  now  act  as  normal  stresses  due  to  the  new 
orientation  of  the  element  caused  by  an  incremental  rotation  during  the  time  increment. 

The  axial  stress  has  the  same  form  as  the  radial  stress,  and  the  tangential  (hoop)  stress  is  also 
similar  except  there  is  no  contribution  from  rotated  shear  stresses. 

The  trial  value  of  the  shear  stress  is  formulated  in  similar  manner. 

=  'tL  +GtxzAt+(a'  -a^CD^At  (26) 

Equations  23—26  assume  an  elastic  response  of  the  material.  If  the  strength  of  the  material  is 
exceeded,  then  plastic  flow  (or  fracture)  will  occur.  The  Von  Mises  yield  criterion  is  used  to 
determine  an  equivalent  stress,  a ,  that  can  be  compared  to  the  uniaxial  tensile  (or  compressive) 
strength  of  the  material.  The  general  form  of  the  equivalent  stress  is 
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Using  deviator  stresses  (instead  of  total  stresses),  Equation  27  can  be  rewritten  as 


If  a  is  not  greater  than  the  equivalent  tensile  strength  of  the  material,  a,  the  final  deviator  and 
shear  stresses  are  as  given  in  Equations  23-26.  If  o  is  greater  than  ct,  then  the  stresses  in 
Equations  23-26  are  multiplied  by  the  factor  (a  /  a) .  When  the  reduced  deviator  and  shear 
stresses  are  put  into  Equation  28,  the  result  is  always  o  =  a .  This  is  known  as  the  radial  return 
algorithm.  The  various  material  strength  models  for  c  are  presented  in  Section  6. 

During  plastic  flow,  it  is  sometimes  necessary  to  determine  the  equivalent  plastic  strain  for  strain 
hardening  effects  on  the  strength  of  the  material  or  to  determine  if  the  material  has  failed.  The 
first  step  in  this  process  is  to  adjust  the  total  strain  rates  to  plastic  strain  rates  by  subtracting  out 
the  elastic  portion  of  the  strain  rates. 


K  =  Cx  -{sT  -  <  +2®,„xl,At)/2GAt 

(29) 

e^  =  e,  -  (sr"‘  -  s*  -  2(o^x^^At)  /  2GAt 

(30) 

eS=ee-(sr-4)/2GAt 

(31) 

YL  =  Yxz - bT  ®xz(sx  - sl)At]/ GAt 

(32) 

Again,  the  other  two  plastic  strain  rates  and  y^)  are  zero  for  this  geometry.  The  general 
expression  for  the  equivalent  plastic  strain  rate  is 
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The  equivalent  plastic  strain,  Ep ,  is  then  obtained  by  integrating  with  respect  to  time. 


pt+At 


=  e*+epAt 


(34) 


After  the  element  stresses  are  obtained,  it  is  necessary  to  determine  concentrated  forces  to  act  on 
the  concentrated  mass  at  the  nodes.  This  is  done  by  obtaining  the  concentrated  forces  which  are 
statically  equivalent  to  the  distributed  stresses  in  the  elements.  The  radial  and  axial  forces  acting 
on  node  i  of  an  element  are 


F*  =  -jux[bia, +CiT,,]--7cAae 


(35) 


F‘=-7tx[cia^  +  bjT^]  (36) 

Note  that  F^  +  Fj  +  F™  =  0  for  each  element,  and  this  ensures  equilibrium  for  the  system.  This 
results  in  conservation  of  momentum  if  there  are  no  external  forces  or  restraints.  The  final  net 
forces  on  node  i  are 


F*=EFi 


(37) 

(38) 


It  is  also  possible  to  add  external  forces  through  applied  pressures. 

The  heat  conduction  algorithm  can  be  obtained  straightforwardly  fi-om  the  ID  heat  conduction 
algorithm  (subsection  3.1.3)  and  the  strain  rate  equations  in  this  subsection.  It  is  reported  in 
Reference  7.  Substituting  temperatures  for  velocities  in  Equation  10  gives 

T  =  ^[(a;  +  b.x + Ciz)Ti  +  (a .  +  b^x  +  c  +  (a„  +  b„x + c„z)t„  ]  (39) 

The  instantaneous  heat  flows  in  the  x  and  z  directions  can  then  be  obtained  from 

q.  =  -k^=  -k(b,T,  +biTj  +b„T.)/2A  (40) 
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3T 

=  -k—  =  -k(cJi  +  c  jT.  +  )  /  2  A 


The  incremental  increase  in  thermal  energy  at  the  nodes  can  be  obtained  by  integrating  the  heat 
flow  with  respect  to  area  and  time. 

AQi  =  tcx(q,bi  +  q^Cj  )At  +  AQ,  (42) 

AQj  =  JCx(q,bj  +q,Cj)At+ajAQ,  (43) 

AQn,  =  Jcx(q,b„  +  q,c„)At  +  a„AQ,  (44) 

Note  that  the  internal  energy  generated  in  the  element  during  the  previous  time  increment,  AQe,  is 
distributed  to  the  nodes  through  aj,  aj,  am  as  defined  by  Equation  2.  The  energy  is  distributed  to 
the  nodes  in  the  same  manner  as  the  mass  is  distributed. 

After  Equations  39-44  have  been  applied  to  all  elements,  the  updated  temperatures  of  the  nodes 
have  the  form 

='!;•+ 1  AQi /MiCpi  (45) 

where  T/'*’^'  and  T/  are  the  temperatures  of  the  nodes  at  times  t+At  and  t,  ZAQj  is  the  sum  of  the 
incremental  heat  contributed  by  all  elements  that  contain  node  i,  M;  is  the  total  mass  of  node  i, 
and  Cpi  is  the  specific  heat  of  node  i. 

The  internal  energy  in  an  element  can  be  used  to  compute  element  pressures.  To  account  for  the 
flow  of  internal  energy  through  the  grid,  the  element  temperature  is  assumed  to  be  the  average  of 
the  nodal  temperatures. 


T  =  (Ti-l-Tj+T„)/3 


The  internal  energy  (per  initial  volume)  is  given  by 


E.  =  (T-T.)p.c, 
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where  To  is  the  initial  temperature  (where  Es  =  0),  po  is  the  initial  density  and  Cp  is  the  specific 
heat  of  the  element  material. 

The  integration  time  increment  must  also  be  boimded  to  ensure  that  the  computations  remain 
stable  for  heat  conduction  (References  6  and  7).  The  heat  conduction  portion  requires 

At<pCph^„/4k  (48) 

where  hmin  is  the  minimum  altitude  of  the  element  and  the  other  terms  recently  have  been 
defined.  This  is  analogous  to  the  time  increment  restriction  for  wave  propagation  in  Equation  4 
of  subsection  3.4.1.  Unless  hmin  is  very  small,  the  wave  propagation  restriction  is  much  more 
severe  than  the  heat  conduction  restriction. 

3.4.5  Quad  Elements 

The  quad  element  algorithm  is  summarized  in  Figure  9.  It  is  simply  an  average  of  the 
components  of  two  quads  composed  of  two  triangles  each. 

The  mass  at  node  i  is  half  the  mass  of  node  i  from  triangular  elements  ijm,  imp,  and  ijp. 

Mi  =  Po[arv;:^"’  +ar'’V'"’‘’  (1) 

where  po  is  the  initial  density  of  the  material  in  the  quad  element,  af”  is  the  fraction  of  the  mass 
of  triangle  ijm  that  is  distributed  to  node  i,  and  is  the  initial  volume  of  triangle  ijm.  Refer  to 

Equations  1-4  in  subsection  3.4.4. 
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Hourglass  Deformation 


Hourglass  Resistance 


Figure  9.  Description  of  the  2D  Quad  Element 
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The  strain  rate  in  the  x  direction  is 


where  A  is  the  area  of  the  quad  element,  zjm  =  zj  -  Zm,  Zmi  =  Zm  —  Zj,  and  U; ,  Uj ,  ,  lip  are  the  x 

velocities  of  the  nodes.  This  result  is  obtained  by  determining  the  x  strain  rates  of  triangles  ijm 
and  imp  (per  Equation  1 1  of  subsection  3.4.4)  and  then  weighting  the  individual  triangle  strain 
rates  by  the  areas  of  the  individual  triangles.  The  same  result  is  obtained  by  considering  the 
alternate  arrangement  of  triangles  ijp  and  jmp. 

Using  this  same  approach,  the  other  strain  rates  and  rotational  rate  can  be  determined. 

ez=^[xjp(v.-Vi)+x„.(vp-Vj)] 

Yxz  =^[xjp(u„  -Ui)+X„i(up  -Uj)  +  Zjp(vi  -  Vj)+Z„i(vj  -Vp)] 

£e  =(ui+Uj+U„+Up)/(x,+Xj+X,„+Xp) 

COxz  =;^[zjp(vi  - v„,)+z^i(zj -Zp)-Xjp(u„ -Ui)-X„i(up -Uj)] 

The  forces  are  determined  in  the  same  manner  as  the  masses.  The  forces  at  node  i  are  half  the 
forces  at  node  i  from  triangles  ijm,  imp,  and  ijp. 

F;  =  -f  +  +Sli,(zjp°.  +Z,jt„)-<J.A]  (7) 

Fz  “  ~'j[zijin(Xn,j<7i  +Zj„X^J+Xi„p(Xp„az  +Zn,pXzz)  Xyp(Xp,<Jp  +  ZjpXpp)]  (8) 

where  Xjj^  is  the  average  x  coordinate  of  triangle  ijm,  and  the  other  terms  are  as  defined 
previously. 


(3) 

(4) 

(5) 

(6) 
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There  is  an  additional  force  required  for  quad  elements  that  is  not  required  for  triangular 
elements.  The  lower  half  of  Figure  9  shows  how  hourglassing  can  occur  in  a  quad  element.  This 
mode  of  deformation  can  occur  without  inducing  any  normal  or  shear  strain  rates.  It  is  detected 
by  evaluating  the  difference  in  the  rotational  rates  of  the  opposite  sides  of  the  quad. 

The  rotational  rates  of  the  four  sides  are 


©mp  =  [(Vm  -  Vp)Xp„  +(up  -  U„)zp^]^p„ 

(9) 

©p.=[(vp-Vi)xpi+(ui-Up)zpi]fpi 

(10) 

(11) 

©mj  =  [(vm  -  Vj)xj„  +(u.  -U„)zj„]£„. 

(12) 

where  is  the  length  of  side  ij  and  the  other  terms  are  as  defined  previously. 

A  viscous  stress  is  used  to  resist  the  hourglassing  (Reference  19)  as  shown  in  Figme  9.  The 
stress  varies  from  ±  qmax  on  each  of  the  opposite  sides.  For  sides  ij  and  mp 

qrr=c„cj3h„,.(<i>.,+(o„)  (13) 

where  Ch  is  an  input  hourglass  viscosity  coefficient,  Cs  is  the  second  velocity  of  the  material,  p  is 
the  density  of  the  material  and  hmin  is  a  characteristic  minimum  length  (area/largest  diagonal). 

For  the  other  two  sides 

=  CHCsph„i„(®pi  +  ®mj)  (14) 

These  stresses  are  then  converted  to  equivalent  nodal  forces  using  the  form 

F',_,=-nx(q»;7x5+qi-'-x,)/3  (16) 
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Both  the  hourglass  forces  and  the  internal  stress  forces  of  Equations  7  and  8  are  then  added  to  the 
appropriate  nodes. 

The  heat  conduction  algorithm  for  quad  elements  is  derived  from  the  triangular  element 
algorithm  in  a  similar  manner.  The  heat  flow  in  the  two  directions  is 

q.=k|^  =  -k[z>,-T„)+z^(Tj-T,)]/2A  (17) 

1.  =  =  -T,)+x^(t,  -T,)]/2A  (18) 

and  the  incremental  increase  in  thermal  energy  at  node  i  is 

AQi=  f  [xij„(zj,q,  +x„jq,)  +  Xi„p(z„pq,  +Xp„q,)+ Xyp(zjpq,  +XpjqJ]At  +  ajAQ,  (19) 

where  Xy„  ,Zj„  ,x„j,  etc.,  are  as  defined  in  Equation  7,  At  is  the  integration  time  increment,  AQe 

is  the  internal  energy  generated  by  the  element  during  the  previous  time  increment,  and  ai  is  the 
fraction  of  AQe  distributed  to  node  i.  The  temperature  and  internal  energy  update  are  similar  to 
those  described  for  the  triangular  elements  in  subsection  3.4.4. 

3.4.6  Bending  Shell  Elements 

A  2D  axisymmetric  bending  shell  element  is  shown  in  Figures  10  and  11.  The  concepts  for  this 
element  are  provided  in  Reference  20.  The  approach  taken  is  to  consider  the  bending  shell  to  be 
represented  by  three  or  five  layered  elements  that  share  the  same  two  nodes.  This  subsection 
provides  the  algorithm  for  the  three  layer  shell,  and  the  five  layer  algorithm  is  similar. 

Bending  moments  about  the  two  end  nodes  are  induced  by  allowing  the  in-plane  stresses  in  the 
outer  layers  to  be  different  from  one  another.  This  is  accomplished  by  allowing  the  in-plane 
strain  rates  in  the  outer  layers  to  be  different  from  one  another.  These  strain  rates  are  dependent 
on  the  nodal  velocities,  as  well  as  the  nodal  rotational  rates  and  the  offset  from  the  center  layer. 

This  algorithm  is  valid  for  small  volumetric  strains  and  thin  shells  (where  the  thickness  is  much 
less  than  the  length).  These  assumptions  lead  to  approximately  equal  lengths  (l,  =  Lj  =  L3  «  l)  , 
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thicknesses  (t,  =  T2  «  T3  «  T  /  3) ,  and  offsets  (h,  =  H2  =  H  =  T  /  3)  for  each  of  the  three  layers. 

Note  that  T  is  the  total  thickness  of  all  three  layers  and  that  Ti,  T2,  and  T3  are  the  thicknesses  of 
the  individual  layers. 

The  shell  element  at  the  top  of  Figure  10  is  in  a  general  orientation  at  an  angle  (j)  from  the  x  axis. 
The  algorithm  is  performed  in  the  x'  -  z'  coordinate  system  as  shown  in  the  lower  portions  of 
Figure  10. 

The  masses  at  nodes  i  and  j  are 

(1) 

Mj  =  ajP„V<,  (2) 


where  po  in  the  initial  density  of  the  element  and  Vo  is  the  initial  volume  of  the  composite 
element  (all  three  layers).  The  fractions  of  the  initial  mass  that  are  assigned  to  nodes  i  and  j  are 


ttj  =(2Xi +Xj)/3(xi +Xj)  (3) 

aj=1.0-aj  (4) 


The  initial  volume  is 

V„  =27txA„  =27CxLJ„ 


(5) 


where  x  is  the  average  initial  x  coordinate  of  nodes  i  and  j,  A©  is  the  initial  area,  Lo  is  the  initial 
length,  and  To  is  the  initial  total  thickness. 

The  mass  moment  of  rotational  inertia  is  required,  and  it  is  expressed  as 

i;  =Mi(L^ +T^)/24  (6) 


Note  that  the  current  (not  initial)  values  of  the  length  and  thickness  (L  and  T)  are  used,  and  this 
requires  I  y  to  be  updated  every  cycle. 
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Figure  10.  Description  of  the  2D  Bending  Shell  Element 
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Stresses 
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Figure  11.  Stresses,  Forces,  and  Moments  for  the  2D  Bending  Shell  Element 
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The  offset  of  the  center  of  layer  2  (above  the  center  of  layer  1),  as  well  as  the  offset  of  layer  3 
below  layer  1,  is  expressed  as 

H  =  V„/67txL  (7) 

where  Vo  is  the  initial  volume  of  the  composite  element,  x  is  the  average  x  coordinate  of  the 
nodes,  and  L  is  the  current  length  of  the  element  (distance  between  nodes).  This  is  a  good 
approximation  if  the  volumetric  strains  are  small  and  the  offsets  are  small  (Axjj  »  ^j2  >  ^i3  > 

AXjj  in  Figure  lo) .  H  is  also  the  thickness  of  each  of  the  layers. 

The  rigid  body  rate  of  rotation  (clockwise  positive)  of  the  composite  element  is 

(j)  =  (z'-Zj)/L  (^) 

where  the  in-plane  (xf  and  x')  and  normal  (z-  and  z')  velocity  components  of  the  center  layer 
are 


x[  =  Zj  sin<t)  +  Xi  cos<|) 

Xj  =  z  j  sin  4)  +  X  j  cos<l) 

z'  =  Zj  cos(t)-Xi  sin(j) 

Zj  =Zj  cos(t)-Xj  sin(|) 

The  nodal  velocities  at  nodes  i  and  j  are  X; ,  Zj 


(9) 

(10) 

(11) 

(12) 

Xj,  and  Zj. 


The  in-plane  velocities  of  the  top  and  bottom  layers  can  now  be  determined. 


x',=x'+H(e,+(i>) 

(13) 

x',=x:-H(e,+4.) 

(14) 

x;,=x'+H(ej+<i>) 

(15) 
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x;3  =  x'-H(0j+<i)) 


(16) 


where  the  in-plane  (x')  velocities  at  node  i  are  x- ,  x'2 ,  and  x'3 ,  the  rotational  velocity  of  node  i 
is  6j ,  and  the  rigid  body  rotational  velocity  of  the  element  is  ^ .  For  nomenclature  consistency, 
the  nodal  x  velocities  will  be  designated  u'jU'j,  and  u'3,  and  the  z  velocities  will  be  designated 
v',v[2,  and  v[3 .  The  velocities  at  node  j  have  similar  notations. 


The  normal  velocities  of  the  three  layers  are  equal  at  node  i  (v'  =  v'2  —  node  j 


The  in-plane  strain  rates  in  the  x'  direction  are 


(17) 

e;2=(u;2-u'2)/L  (18) 

ex3  =  (uj3-u;3)/L  (19) 

The  hoop  strain  rates  are 

£9,  =h, /x  (20) 

e92=U2/x  (21) 

^63  =^3/^  (22) 


where  u, ,  U2 ,  U3  are  the  average  x  velocities  of  the  two  ends  of  the  three  shell  layers,  and  x  is 
the  average  x  coordinate  of  the  composite  element. 


The  shear  strain  rates  for  all  three  layers  are  identical 


'dV  .  ^ 

,ax  V 


V'-V'  0i-he: 
— — + - - 


(23) 
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where  v'  =  z'  and  v'  =  z'  are  from  Equations  11  and  12. 


After  the  in-plane  strain  rate,  e' ,  the  hoop  strain  rate  ,  and  the  shear  strain  rate,  ,  have 
been  determined,  the  through-thickness  strain  rate,  e,,  must  be  determined.  The  plane  stress 
condition  of  the  shell  requires  the  through-thickness  stress  to  be  CT^  =  0 .  Therefore,  it  is 
necessary  to  iterate  to  find  the  strain  rate  that  produces  the  correct  stress  (a'  =  o) . 

The  iteration  algorithm  is  shown  in  Figure  12.  The  initial  through-thickness  strain  rate  and  stress 
[e'(l)  and  a'(l)]  come  from  assuming  incompressible  flow  during  the  previous  cycle.  The  two 

in-plane  normal  strain  rates  and  89)  are  obtained  in  the  standard  manner,  and  then  the 

through-thickness  strain  rate  is  set  to 

e;(l)  =  -(£;  +  £») 


which  results  in  incompressible  flow  for  the  cycle  {e'  +  e'  +  e,  —  o) . 

From  the  current  strain  rates,  the  previous  stresses,  and  the  previous  volumetric  strain,  the  net 
through-thickness  stress  is 

<  =  s;-P  (27) 

where  s'  is  the  deviator  stress,  determined  from  the  current  strain  rates  and  the  previous  stresses. 
The  pressure  is  determined  from  P  =  -Key  where  K  is  the  bulk  modulus  and  Sy  is  the  previous 
volumetric  strain.  The  strain  rate  and  stress  for  the  incompressible  condition  are  represented  by 
e'(l)  and  a'(l)  in  Figure  12. 
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Figure  12.  Normal  Stress  Iteration  Algorithm  for  2D  and  3D  Bending  Shells,  and 
Plane  Stress  Geometry 

The  next  step  is  to  determine  e'(2)  and  cf'(2)  ,  as  shown  in  Figure  12.  This  condition  is 
assumed  to  be  that  which  would  result  from  a  compressible  elastic  strain  rate  of 

E;(2)  =  -v(e;  +  E,)/(l-v)  (28) 

where  v  is  Poissons  ratio,  and  e'  and  are  identical  to  those  of  Equation  26. 

Using  the  new  strain  rate  it  is  possible  to  determine  a  new  normal  stress  from  Equation  27.  Both 
the  deviator  stress  and  the  pressure  are  different  from  the  incompressible  step.  The  pressure  is 
different  because  the  volumetric  strain  has  been  modified  to: 

£.(2)  =  £.(l)+[£;(2)-£;(l)]At  (29) 

where  At  is  the  integration  time  increment. 


A  t19616.doc 


57 


Now  by  linear  interpolation  it  is  possible  to  determine  £^(3)  for  —  0 .  The  new  computed 
value  of  a; (3)  will  generally  be  close  to  =  0,  but  an  additional  iteration  may  be  required. 
This  is  represented  by  e'(4) ,  and  it  is  the  average  e'  computed  from  interpolation/extrapolation 

of  points  1  and  3,  and  2  and  3. 

The  resulting  stresses  are  shown  on  the  upper  half  of  Figure  11.  These  stresses  are  integrated 
over  the  appropriate  areas  to  obtain  the  forces  on  the  lower  half  of  Figure  11. 

For  layer  1,  the  in-plane  force,  the  shear  force,  and  the  bending  moment  are  given  by 


f;=f;=2jcxh<, 

(30) 

v;=v;=27txHx::, 

(31) 

M‘  =-V,L/2 

(32) 

Mj  =  V,L/2 

(33) 

where  x  is  the  average  x  coordinate  of  nodes  i  and  j,  H  is  the  thickness  of  the  individual  element 
layer  from  Equation  7,  and  L  is  the  current  length.  These  forces  do  not  include  the  force  due  to 
the  hoop  stress.  The  in-plane  forces,  the  shear  forces,  and  the  bending  moments  for  layers  2  and 
3  are  determined  in  a  similar  manner.  The  moment  component  (mJ  and  M])  is  included  to  keep 

the  element  in  equilibrium. 

The  net  forces  and  moments  on  nodes  i  and  j  (from  all  three  layers)  are  as  follows: 


F;.  =  Fi  =  F,+F,+Fj 

(34) 

F;.=Fj.=(v,+V3  +  V,) 

(35) 

M‘ =  M‘ +  M’ +  M‘ +  FjH  -  FjH 

(36) 

M^  =  Mj+M^  +M^  -t-FjH-FjH 

12  3-^  J 

(37) 
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These  are  per  the  sign  convention  at  the  top  of  Figure  1 1,  and  in  the  element  (local)  coordinate 
system. 

The  net  forces  on  node  i,  in  the  system  coordinate  system,  are 

=  -F^.  cos<l)  +  F^,  sin(j)-7tA(a0,  +a02  +093)  (38) 

Fz  =-Fx'Sin(|)-F^,cos<t)  (39) 

where  A  is  the  cross-sectional  area  of  each  individual  shell  layer  and  Gg, ,  Gqj  ,  and  0^3  are  the 

hoop  stresses  for  each  of  the  three  layers.  The  net  nodal  force  on  node)  is  determined  in  a 
similar  manner. 

An  assessment  of  the  accuracy  of  the  three  and  five  layer  shell  algorithms  is  shown  in  Figure  13. 
The  errors  are  1 1  percent  for  the  three  layer  algorithm  and  4  percent  for  the  five  layer  algorithm. 

The  heat  conduction  algorithm  is  similar  to  the  ID  algorithms,  except  the  cross-sectional  area  (in 
the  x'  direction)  is 


A  =  27CxH  (40) 

where  X  is  the  average  x  coordinate  of  the  nodes  and  H  is  the  thickness  of  the  individual  layer. 
For  the  heat  conduction  option  all  three  layers  will  have  equal  temperatures  because  all  layers  are 
attached  to  the  same  nodes. 

3.4.7  Membrane  Shell  Elements 

The  2D  membrane  shell  elements  are  a  simplified  form  of  the  bending  shell  elements  presented 
in  subsection  3.4.6  and  shown  in  Figure  10.  The  membrane  shell  element  has  only  in-plane 
stresses  and  is  composed  of  only  a  single  layer.  It  cannot  produce  bending  or  transverse  shear 
stresses.  Furthermore  it  is  assumed  to  be  incompressible,  and  this  eliminates  the  need  to  iterate 
to  obtain  the  through-thickness  stress.  This  algorithm  is  presented  in  Reference  21  for  plane 
stress  geometry. 
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M  =  2/9  aH2  =  .222  oH^ 
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M  =  6/25  aH2  =  .240  oH^ 


4%  Low 


Figure  13.  Accuracy  Assessment  for  Bending  Sheiis 

Referring  back  to  Equation  23  of  subsection  3.4.6,  the  through-thickness  strain  rate  for 
incompressible  deformation  is 

ez  =  -(ex+ee) 

where  e'  is  the  in-plane  strain  rate  and  is  the  hoop  strain  rate. 

The  deviator  stresses  (s' ,s' .s^)  can  then  be  determined  from  e'  ,  e' ,  £9 ,  which  give  the  net 
stresses 
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a:=s:-p 


a:=s:-p 


a;=s;-p 


(2) 

(3) 

(4) 


In  Equation  3,  setting  a'  =  0  gives  P  =  s' ,  and  this  provides  an  explicit  determination  of  the 
stresses  in  Equations  2-4. 

3.4.8  Nonreflective  Boundary  Elements 

Nonreflective  boundary  elements  (References  8  and  9)  have  been  available  historically  to  provide 
nonreflective  damping  for  the  modeling  of  large  infinite  bodies  such  as  soil  or  water.  They  are 
intended  primarily  to  absorb  elastic  waves,  but  they  are  also  effective  for  other  applications 
(Reference  10).  The  2D  nonreflective  boundary  elements  are  incorporated  as  two-node, 
infinitely-thin,  massless  elements,  as  shown  in  Figure  14. 
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Figure  14.  Description  of  the  2D  Nonreflective  Boundary  Element 
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The  average  velocity  normal  to  the  boimdary  element  is 


=  u  sin  0  -  V  cos  0 

(1) 

where  h  =  (uj  +  Uj) /  2  is  the  average  velocity  in  the  x  direction,  v  is 
z  direction,  and  0  is  the  angle  between  the  element  and  the  x  axis. 

the  average  velocity  in  the 

The  average  velocity  parallel  to  the  element  (from  node  i  to  node  j)  is 

Vp  =  hcos0  + vsin0 

(2) 

The  total  normal  and  parallel  forces  are 

(3) 

Fp  “Po^shear-^ 

(4) 

where  p„  is  the  initial  density  of  the  material,  A  =  27tx  is  the  surface  area  of  the  element,  and  Vn 
and  Vp  are  the  velocities.  The  longitudinal  and  shear  velocities  are 

C3=^(k,+4G/3)/p„  (5) 

^shear  ~  ^  Po 

(6) 

where  K,  is  the  bulk  modulus  and  G  is  the  shear  modulus. 

Now  the  total  forces  must  be  aligned  with  the  two  principal  axes  and  distributed  to  the  two  nodes 

F-  =  Fj  =  (Fn  sin0  +  Fp  cos0)/2 

(7) 

F’  =  Fj  =  (Fn  cos0+Fp  sin0)/2 

(8) 
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3.5  2D  AXISYMMETRIC  GEOMETRY  WITH  SPIN 


The  2D  axisymmetric  geometry  with  spin  is  identical  to  the  2D  axisymmetry  geometry  described 
previously  in  subsection  3.4,  except  that  this  geometry  allows  the  nodes  to  experience  rotations 
(or  spin)  about  the  axis  of  revolution,  as  shown  in  Figure  15. 


Z 


X 

Figure  15.  Description  of  the  2D  Axisymmetric  Triangular  Eiement  with  Spin 

For  this  geometry  the  concentrated  masses  at  the  nodes  can  be  visualized  as  concentric  circular 
rings  contained  in  planes  that  are  perpendicular  to  the  axis  of  revolution.  These  rings  can  move 
up  and  down  along  the  axial  direction,  they  can  expand  and  contract  in  the  radial  direction,  and 
they  can  rotate  about  the  axis  of  revolution.  This  additional  spinning  degree  of  freedom 
introduces  centrifugal  forces,  and  also  allows  for  two  additional  shear  strain  rates  and  shear 
stresses.  The  algorithm  is  reported  in  References  3  and  1 1 . 
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The  three  normal  strain  rates  (e^,e^,ee)  and  the  shear  strain  rate  in  the  x-z  plane  (y„)  for  a 
triangular  element  are  identical  to  those  described  in  subsection  3.4.  The  two  addition  shear 
strain  rates  are 


Yxe  = 


dw  w 
dx  X 


(1) 


Yze  = 


dw 

dz 


(2) 


where  Wj  =  XjBj  is  the  tangential  velocity  of  node  i  under  a  rotational  rate  0; .  The  derivatives 
^  3  w  3  w 

—and—  are  obtained  in  the  same  manner  as  the  other  derivatives  were  obtained  in 
ydx  dz  J 

subsection  3.4.  In  Equation  1,  W  is  the  average  tangential  velocity  of  the  three  nodes  and  x  is 
the  average  x  coordinate. 

The  shear  strain  rates  in  Equations  1  and  2  lead  to  two  new  shear  stresses  (x^ean^T^e),  that  are 
obtained  in  the  same  manner  as  described  in  subsection  3.4.  The  resulting  tangential  force  is 


Fe  =  -7tx| 


(3) 


where  bj  =  zj  -  Zm  and  Ci  =  Xm  -  xj,  as  defined  previously.  The  factor 
equilibrium  in  the  9  direction. 


is  required  to  maintain 


The  equations  of  motion  are  modified  and  expanded  for  this  geometry.  The  radial  acceleration 
(given  previously  in  subsection  3.4.1)  is  expanded  to  include  centrifugal  force  due  to  spin. 

x‘ =  FJ /Mi+x‘(e‘"f  (4) 

where  Fj  and  M;  are  the  force  and  mass  described  in  subsection  3.4.1,  x\  is  the  radial 
coordinate  at  time  =  t  and  0*‘  is  the  rotation  velocity  before  it  is  updated  at  time  =  t. 
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The  updated  radial  velocity  is 


x*"^  =  (x|  +  X J  At)  (l  -  Cp  At  /  xj )  (5) 

where  x|”  is  the  constant  velocity  for  the  previous  time  increment  and  At  is  the  average  of  the 
two  integration  time  increments  about  time  =  t.  The  expression  in  the  second  set  of  parentheses 
can  be  used  to  damp  out  the  radial  velocities  to  give  steady-state  solutions  for  spinning  bodies.  If 
the  constant,  Cd,  is  set  equal  to  twice  the  soimd  velocity  of  the  material,  the  system  will  be 
approximately  critically  damped  and  the  steady-state  solution  will  be  rapidly  attained. 

For  the  rotational  equations  of  motion  it  is  necessary  to  consider  the  angular  momentum,  H,  of 
each  concentrated  mass.  This  gives 

=H*"+x‘ f;  At  (6) 

By  substituting  Hj  =  0j  xf  Mj  into  Equation  (6),  it  is  possible  to  determine  the  updated  rotation 
velocity. 


^  xJ  At 

1  .  * 

r  F*  ^ 

Te 

(xr)^ 

ImJ 

(7) 


It  should  be  noted  that  even  if  the  net  circumferential  force,  F^  is  equal  to  zero,  it  is  possible  for 

the  spin  to  change  if  the  radius  changes  between  times  t  and  t+At.  It  is  therefore  necessary  to 
obtain  the  new  radial  position  at  t+At  before  obtaining  the  new  rotation  velocity  at  t+At. 

3.6  2D  PLANE  STRAIN  GEOMETRY 

The  2D  plane  strain  geometry  is  a  simplified  form  of  the  2D  axisymmetric  geometry.  The 
primary  difference  is  that  the  plane  strain  geometry  is  for  a  unit  thickness  and  has  no  hoop  strains 
or  strain  rates. 

Referring  back  to  subsection  3.4.4  and  Figure  4  for  2D  axisymmetric  triangular  elements,  some 
of  the  revised  equations  are  as  follows: 
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The  mass  at  node  i  is 


3.7  2D  PLANE  STRESS  GEOMETRY 

The  2D  plane  stress  geometry  is  similar  to  the  2D  plane  strain  geometry.  The  difference  is  that 
for  plane  strain  geometry  all  of  the  strain  is  in  the  x-z  plane  and  there  is  no  strain  normal  to  the 
plane  =  Ey  =  Oj .  For  plane  stress  geometry,  all  the  stress  is  in  the  x-z  plane  and  there  is  no 
stress  normal  to  the  plane  (Cy  s  Oj .  The  initial  plane  stress  geometry  has  an  initial  thickness  of 
Ayo  =  1.0,  but  the  thickness  can  change  as  the  grid  is  deformed.  The  original  algorithm  was  for 
incompressible  flow  (Reference  21),  but  the  current  algorithm  includes  compressibility. 
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The  plane  stress  algorithm  is  identical  to  the  plane  strain  algorithm,  except  for  the  determination 
of  the  through-thickness  strain  rate,  and  the  determination  of  the  forces.  Because  the  through¬ 
thickness  stress  must  be  Gy  =  0  it  is  necessary  to  iterate  to  determine  this  stress.  This  is 
accomplished  by  using  the  iteration  algorithm  described  in  subsection  3.4.6  and  shown  in 
Figure  12.  For  bending  shells  the  known  strain  rates  are  e'  and  ,  and  the  through-thickness 
strain  rate  (to  be  determined)  is  e' .  For  plane  stress  geometry,  the  known  strain  rates  are  and 
,  and  the  through-thickness  strain  rate  (to  be  determined)  is  .  Other  than  the  change  in 
nomenclature,  the  algorithms  for  the  two  cases  are  identical. 

The  forces  are  slightly  modified  from  those  of  the  plane  strain  geometry 


Fx=-Ay(bia,-hCiX^)/2 

(1) 

Fz  =  -Ay(cia,-i-bjX„)/2 

(2) 

where  Ay  is  the  current  thickness  and  the  other  terms  are  identical  to  those  of  the  plane  strain 
geometry  given  in  Equations  4  and  5  of  subsection  3.6.  The  current  thickness  is 

Ay  =  V/A  =  V„(l-HeJ/A  (3) 

where  V  is  the  ciurent  volume,  A  is  the  current  area,  Vo  is  the  initial  volume,  and  Cy  is  the  current 
volmnetric  strain. 

3.8  3D  GEOMETRY 

A  description  of  3D  geometry  is  shown  in  Figure  16.  This  figure  shows  a  tetrahedral  element, 
but  brick  elements,  bar  elements,  bending  shell  elements,  membrane  shell  elements,  and 
nonreflective  boimdary  elements  are  also  available.  The  early  developments  of  this  work  are 
reported  in  References  2, 4, 14,  and  22-24. 
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T19616 

Figure  16.  Description  of  the  3D  Tetrahedral  Element 
3.8.1  Equations  of  Motion 

The  accelerations,  velocities,  and  positions  are  determined  as  follows. 
x'=F,‘/Mj 
y*=F;/Mi 
z‘=Fj/Mi 
x*-"  =x*'+x*At 
yr=yr+y.‘^t 
=zr+z‘^ 
xr^‘  =x‘+x|^At 
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(8) 


yr=y!+yrAt 

z*^^' =z|+z‘^At  (9) 

Equations  1-3  provide  the  x,  y,  and  z  accelerations  at  time  =  t.  F^' ,  Fy ,  and  F^'  are  the  net  x,  y, 

and  z  direction  forces  on  node  i.  These  forces  are  contributed  to  by  all  of  the  elements  that 
contain  node  i,  and  are  obtained  from  the  element  computations  in  the  previous  cycle  of 
integration.  Mj  is  the  total  mass  at  node  i,  and  it  contains  a  fraction  of  the  mass  of  all  elements 

that  contain  node  i.  Although  Figure  16  shows  only  a  tetrahedral  element,  a  node  can  have 
various  element  types  attached  to  it.  This  means  it  is  possible  for  an  individual  node  to  have 
masses  and  forces  from  different  element  types.  The  specific  algorithms  for  masses  and  forces 
are  provided  in  subsections  3. 8.4-3. 8. 9. 

The  updated  velocities  in  Equations  4-6  are  constant  for  the  interval  between  time  =  t  and  time  = 
t+At.  The  constant  velocities  for  the  previous  time  increment  are  x*” ,  y •” ,  and  zj“ ,  and  At  is 

the  average  of  the  two  integration  time  increments  about  time  =  t. 

The  integration  time  increment  is  limited  to 

At  =  C,[h/(V?+7g^  +c^ )]  (10) 


where  g^  =  CqQ/p,  h  is  a  characteristic  length  (described  later  for  each  element  type),  and  Cs  is 
the  sotmd  velocity  of  the  material  (References  2  and  5).  Q  is  the  artificial  viscosity,  and  Cq  is  a 
constant  for  the  artificial  viscosity.  These  are  described  in  Section  6.  The  Courant  soimd  speed 
fraction,  Q,  must  be  less  than  unity  (Q  <  1 .0)  to  ensure  numerical  stability.  Q  =  0.9  is  a  typical 
value. 

For  3D  bending  shell  elements  there  are  also  rotational  degrees  of  freedom  on  each  node.  The 
rotational  velocities  of  node  i  are  updated  in  the  same  general  manner  as  the  translational 
velocities  in  Equations  1-6,  except  that  they  are  coupled,  and  therefore  more  complex 
(Reference  25).  The  system  of  equations  to  be  solved  is  as  follows: 


B  t19616.doc 


69 


(11) 


"yx 


^xy  ^xz 
lyy  lyz 

Izy  Izz 


M.  +(i„  -  i„Ke. + ('.>».  -  -  eJ) 

M, +(u  -i„)e,e. +(i,i  -i^e,)e, +i„(6j  -eO 

M.  +  (l„-I„)eA+My-‘AK  +  >.y(95-80. 


The  rotational  accelerations  0y,  ,  the  moments  of  inertia  (Ixx  •••Izz)  and  the  net  moments 

(m^,  My,  are  acting  on  node  i  at  time  =  t.  The  rotational  velocities  6y,  are  for  the 
time  increment  prior  to  time  =  t. 

The  updated  accelerations  are  determined  from  Equation  1 1  and  are  designated  0^,  0y ,  0jj .  The 
updated  velocities  for  the  next  time  increment  are 


0‘;=0‘-+0‘At 


0;"=0;-+0;At 


0'"=0‘-+0'At 


(12) 

(13) 

(14) 


The  rotational  displacements  (0^,  0y,  0^  are  not  required. 

3.8.2  Sliding  interfaces 

The  3D  sliding  interface  algorithm  is  similar  to  the  2D  sliding  algorithm  described  in  subsection 
3.4.2,  but  it  does  not  have  the  order  independent  algorithm  for  automatic  sliding.  Although  an 
automatic  sliding  algorithm  is  available  in  3D,  the  searching  time  can  be  excessive  and  it  is 
usually  more  computationally  efficient  for  the  user  to  specify  the  master  and  slave  regions. 
Although  the  3D  sliding  algorithm  has  been  improved  in  recent  years,  the  early  algorithms  are 
reported  in  References  14  and  22-24. 

The  interface  determination  algorithm  and  the  search  algorithm  will  not  be  described  herein, 
except  to  note  that  they  are  analogous  to  the  2D  algorithms  described  in  subsection  3.4.2. 

Figure  17  shows  a  slave  node  in  contact  with  a  triangular  master  surface  defined  by  nodes  i,  j, 
and  k.  The  normal  velocities  of  the  slave  node  and  the  three  master  nodes  are  adjusted  first. 
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Three  of  the  foxir  conditions  required  to  determine  these  new  velocities  involve  conservation  of 
linear  momentum  normal  to  the  master  plane  and  conservation  of  angular  momenta  about  the  two 
axes  in  the  plane.  These  conditions  define  the  fractions  (Rj,  Rj,  Rk)  of  linear  momenta  transferred 
from  the  slave  node  to  the  three  master  nodes.  The  resulting  momenta-conserving  relationships 
between  the  instantaneous  normal  velocity  changes  of  the  slave  node  (AV^^  )  and  the  master 

nodes  ,  AVl* ,  AV^ )  have  the  form 


AVj^  =-RiM,AV3^  /Mj 


(1) 


where  Ms  and  Mi  represent  the  mass  of  the  slave  node  and  the  master  node  i. 


Figure  17.  Description  of  the  3D  Siiding  Interface  Contact  Algorithm 
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The  final  condition  involves  equating  the  normal  velocities  of  the  slave  node  and  the  master 
surface  at  the  slave  node  position.  For  this  case,  it  is  necessary  to  determine  the  velocity  of  the 
master  surface  at  the  slave  node  location.  Here  the  velocity  is  determined  within  a  two- 
dimensional  triangle.  The  resulting  master  surface  velocity  in  the  x  direction,  at  the  slave  node 
position,  is 

Um  =qiUi+qA+qkUk  ^ 


The  other  velocity  components  (v„  and  w^ )  have  similar  form  to  Equation  2  and  the  geometry 
constants  (qi,  qj,  qk),  have  the  form 


qi=^hyk-Xkyj+(yj 


(3) 


where  Axy  is  the  projected  cross-sectional  area  of  the  triangular  master  surface  on  the  x-y  plane, 
and  the  other  terms  represent  the  coordinates  of  the  master  and  slave  nodes.  The  projections 
could  also  be  on  the  x-z  or  y-z  plane  depending  on  the  orientation  of  the  master  plane. 


By  equating  the  normal  velocities  and  substituting  the  relationships  of  Equation  1,  the  normal 
velocity  change  imposed  on  the  slave  node  is  defined  as 


Avr  = 


yN  _  yN 


1  +  qjRjM,  /  Mj  +  qjRjM,  /  Mj  +  q^RkM,  / 


(4) 


where  the  slave  node  velocity  normal  to  the  master  surface  is 

=  AUj  +  BVj  +  CWj  (3) 

where  A,  B,  C  are  direction  cosines  of  the  vector  normal  to  the  master  surface.  The  normal 
velocity  of  the  master  surface,  ,  has  the  same  form  as  Equation  5  and  the  velocity  changes  to 

the  three  master  nodes  can  be  defined  from  Equation  1.  The  specific  velocity  changes  in  the  x,  y, 
and  z  directions  are  obtained  by  multiplying  AV  by  the  appropriate  direction  cosines. 
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Frictional  effects  are  constrained  to  the  plane  of  the  master  surface  eind  the  direction  opposes  the 
relative  motion.  The  net  magnitude  of  the  frictional  velocity  change  to  the  slave  node  is 
proportional  to  the  normal  velocity  change  to  the  slave  node,  or 

AVf  =  f,AV,^  (6) 


where  fs  is  the  coefficient  of  friction. 

The  velocity  components  in  the  plane  of  the  master  smface  are  obtained  by  subtracting  the 
velocity  components  normal  to  the  master  plane,  from  the  total  velocity  components.  These 
velocity  components  for  the  slave  node  are  expressed  as 


uf  =u  -AV.^ 


(7) 


v?=v  -BV!' 


(8) 


wf  =w  -CV.^ 


(9) 


The  in-plane  velocities  of  the  master  surface  (u^ »  )  are  obtained  in  a  similar  manner. 

The  relative  velocity  components  (u^, ,  v^, ,  w^, )  are  then  obtained  in  the  form  of 


•  P  •  P  •  P 

u  ,  =  u  —  u 

'^rel  '*in  “s 


(10) 


Finally,  the  friction  induced  velocity  changes  (AUj  ,Av5  ,Aw3 )  to  the  slave  node  have  the  form 


Aii*^  = 

S 


.\P 


u; 


rel 


+(vL)'  +(Ki)\ 


AVJ 


(11) 


The  momentum  change  of  the  master  nodes,  due  to  the  frictional  force,  gives  in-plane  velocity 
changes  to  master  node  i(Auf,  Avf,  Awf)  of  the  form 

Aur  =-RiM,Au3^/Mi  (12) 


B  t19616.doc 


73 


For  the  standard  (not  automatic)  sliding  algorithm  the  slave  nodes  carmot  also  be  master  nodes 
(for  a  given  interface).  Here  the  velocity  changes  are  collected  after  all  slave  nodes  have  been 
processed  and  there  is  no  order  dependence. 

For  the  automatic  sliding  algorithm  the  slave  and  master  node  velocities  are  updated  as  they  are 
processed,  and  this  results  in  some  order  dependence. 

Future  plans  are  to  incorporate  a  new  3D  algorithm  that  is  analogous  to  the  current  2D  algorithm. 

3.8.3  Erosion 

The  three-dimensional  algorithm  that  follows  should  be  applied  only  to  problems  in  which 
erosion  is  the  primary  mode  of  penetration.  It  is  reported  in  Reference  16.  It  is  only  available  for 
tetrahedral  elements.  For  many  problems,  the  master  surface  is  initially  on  the  top  (or  any  other 
side)  of  a  plate  and  it  cannot  involve  other  sides  (except  as  special  cases)  during  erosion.  The 
precise  restriction  for  the  general  algorithm  to  follow  is  that  no  element  with  a  free  surface  can  be 
eroded  unless  that  free  surface  was  initially  designated  as  a  master  surface.  The  two  exceptions 
are  for  a  free  surface  on  the  y  =  0  plane  of  symmetry,  and  for  a  free  surface  on  the  bottom  of  the 
plate. 

Figure  1 8  shows  a  partially  eroded  master  surface  defined  as  a  set  of  adjacent  triangles.  The 
initial  master  surface  included  the  entire  top  surface  of  the  plate.  The  slave  nodes  are  essentially 
above  the  master  surface  and  must  always  view  the  associated  master  triangles  so  that  the  nodes 
of  the  triangles  (Ml,  M2,  M3)  are  in  a  coimterclockwise  order.  The  slave  nodes  can  be  randomly 
ordered  and  are  not  allowed  to  penetrate  through  the  master  surface.  The  master  surface  is 
allowed  to  erode  by  totally  failing  tetrahedral  elements  which  have  one,  two,  three  or  four  sides 
on  the  master  surface.  Totally  failed  elements  do  not  develop  stresses  or  pressures;  they 
essentially  disappear  except  that  mass  is  retained  at  the  nodes. 
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Figure  18.  Description  of  the  3D  Erosion  Aigorithm 

Before  the  computations  begin,  the  preprocessor  must  store  the  groups  of  three  nodes  defining 
the  individual  triangular  master  planes.  It  must  also  permanently  mark  each  node  contained  on 
the  master  surface,  and  it  must  count  and  store  the  number  of  elements  attached  to  each  node.  If 
three  or  four  nodes  of  an  element  are  marked,  then  the  element  is  on  the  master  surface  and  is 
eligible  for  failure  if  it  exceeds  the  specified  erosion  strain.  When  an  element  is  totally  failed,  the 
element  count  for  each  of  the  four  nodes  is  decremented,  so  that  the  cotmt  reflects  the  number  of 
imfailed  elements  attached  to  the  node.  When  the  count  reaches  zero,  the  node  is  made  a  slave 
node  because  it  is  no  longer  associated  with  the  master  sixrface.  This  new  slave  node  (which 
continues  to  have  mass,  momentum,  and  kinetic  energy)  can  then  interact  with  the  master  surface 
but  cannot  pass  through  it. 

The  six  elements  in  Figure  18  (elements.  A,  B,  C,  D,  E,  F)  have  three  or  four  nodes  on  the  master 
surface  and  have  one,  two,  three  or  four  sides  on  the  master  surface.  The  following  describes  the 
changes  necessary  to  redefine  the  master  surface  and  the  slave  nodes: 

.  Element  A  has  one  side  (M1-M2-M3)  on  the  master  surface.  All  four  nodes  (Ml,  M2, 
M3,  and  M4)  are  marked.  Because  only  one  side  (M1-M2-M3)  is  contained  in  the  list 
of  master  triangles,  it  is  known  that  the  other  three  sides  are  not  on  the  master  surface. 

If  the  equivalent  plastic  strain  or  volumetric  strain  of  element  A  exceeds  the  erosion 
strain,  the  element  is  totally  failed.  The  updated  master  siirface  has  one  triangular  side 
(M1-M2-M3)  removed  and  three  triangular  sides  of  adjacent  elements  (M1-M2-M4, 
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M1-M4-M3,  M2-M3-M4)  added.  The  nodes  of  the  new  sides  must  be  given  in  the 
proper  order  so  that  the  slave  nodes  see  a  counterclockwise  ordering.  The  element 
counts  of  nodes  Ml,  M2,  M3,  M4  are  decremented,  but  all  four  nodes  remain  on  the 
master  surface  because  none  of  the  element  coimts  goes  to  zero. 

.  Element  B  has  two  sides  (M5-M6-M7,  M6-M8-M7)  on  the  master  surface.  The  four 
nodes  M5,  M6,  M7,  and  M8  are  marked.  If  the  equivalent  plastic  strain  or  volumetric 
strain  of  element  B  exceeds  the  erosion  strain,  the  element  is  totally  failed.  The 
previous  master  surface  has  two  sides  (M5-M6-M7,  M6-M8-M7)  removed  and  two 
other  surfaces  (M5-M6-M8,  M5-M8-M7)  added.  The  element  counts  of  nodes  M5, 
M6,  M7  and  M8  are  decremented,  but  none  of  the  counts  goes  to  zero. 


.  Element  C  has  three  sides  (M9-M10-M1 1,  M9-M1 1-M12,  M9-M12-M10)  on  the 
master  surface.  The  four  nodes  M9,  MIO,  Ml  1,  and  M12  are  marked.  If  the  equivalent 
plastic  strain  or  volumetric  strain  of  element  C  exceeds  the  erosion  strain,  the  element  is 
totally  failed.  The  previous  master  surface  has  three  sides  (M9-M10-M1 1,  M9-M1 1- 
M12,  M9M12-M10)  removed  and  one  other  surface  (MlO-Ml  1-M12)  added.  The 
element  counts  of  nodes  M9,  MIO,  Ml  1,  and  M12  are  decremented,  but  none  of  the 
coimts  goes  to  zero. 


•  Element  D  has  four  sides  (M13-M14-M15,  M13-M16-M14,  M14-M16-M15,  M13- 
M15-M16)  on  the  master  surface.  The  four  nodes  M13,  M14,  M15,  and  M16  are 
marked.  If  the  equivalent  plastic  strain  or  volimietric  strain  of  element  D  exceeds  the 
erosion  strain,  the  element  is  totally  failed.  The  previous  master  surface  has  all  four 
sides  removed  and  none  added.  The  element  count  of  nodes  M13,  M14,  M15,  and  M16 
are  decremented.  The  element  counts  of  nodes  M13,  M15,  and  M16  do  not  go  to  zero 
and  therefore  remain  on  the  master  surface.  The  element  count  on  node  M14  does  go  to 
zero,  however,  so  it  is  designated  to  be  a  slave  node. 

•  Element  E  is  a  special  case  because  it  has  one  side  (Ml 8,  Ml 9,  M20)  on  a  plane  of 
symmetry  at  y  —  0,  which  was  not  initially  designated  as  a  master  surface.  The  general 
algorithm  would  indicate  one  side  (M17-M18-M19)  to  be  deleted  and  three  sides  to  be 
added  (M17-M18-M20,  M17-M20,  M19,  M18-M19-M20).  The  side  on  the  plane  of 
symmetry  at  y  =  0  (M18-M19-M20)  should  not  be  added,  however.  This  side  is  readily 
identified  and  deleted  by  noting  that  y  =  0  for  nodes  M18,  M19,  and  M20.  For  this  case 
the  element  counts  of  nodes  M17,  Ml  8,  M19,  and  M20  are  again  decremented.  Node 
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M20  was  not  previously  marked  and  must  therefore  be  marked  to  identify  it  as  a  part  of 
the  updated  master  smface.  It  should  be  noted  that  the  general  algorithm  would  handle 
this  special  case  if  the  triangular  element  sides  on  the  plane  of  symmetry  would  initially 
be  designated  as  part  of  the  master  surface.  Because  this  would  significantly  increase 
the  size  of  the  master  surface  (and  the  associated  searching  time),  it  is  instead  treated  as 
a  special  case. 

•  Element  F  is  another  special  case  because  it  is  eroded  through  a  free  surface  (bottom  of 
the  plate)  which  was  not  initially  designated  as  a  master  surface.  Here  the  general 
algorithm  would  remove  two  sides  (M21-M22-M23,  M21-M23-M24)  and  add  two  new 
sides  (M21-M22-M24,  M22-M23-M24).  The  side  on  the  bottom  of  the  plate  (M22- 
M23-M24)  should  not  be  added,  however.  This  side  can  be  readily  identified  if  the 
nodal  numbering  system  proceeds  downward  by  layers,  thus  placing  the  highest  node 
numbers  on  the  bottom  surface  of  the  plate.  Under  these  conditions,  node  numbers 
M22,  M23,  and  M24  would  all  be  greater  than  (or  equal  to)  the  lowest  node  number  on 
the  bottom  of  the  plate,  and  that  triangle  would  therefore  not  be  added  to  the  list  of 
master  triangles.  Other  approaches  could  be  used  for  different  nodal  numbering 
systems.  Again  it  should  be  noted  that  the  general  algorithm  would  handle  this  special 
case  if  the  element  sides  on  the  bottom  of  the  plate  would  initially  be  designated  as  part 
of  the  master  surface. 

To  ensure  that  there  is  no  cross-over  of  material  between  the  two  surfaces  (such  as  a  projectile 
and  a  target),  a  two-step  approach  (two  separate  sliding  interfaces)  is  used  with  the  standard  (not 
automatic)  sliding  algorithm.  For  the  first  interface,  all  of  the  projectile  nodes  and  the  eroded 
target  nodes  are  slave  to  the  master  surface  in  the  target.  For  the  second  interface,  all  of  the 
nodes  in  the  center  of  the  target,  and  the  eroded  projectile  nodes,  are  slave  to  the  master  surface 
on  the  projectile.  The  eroded  nodes  are  therefore  restrained  from  passing  through  either  the 
projectile  or  the  target.  The  automatic  sliding  algorithm  handles  both  surfaces  automatically. 

3.8.4  Tetrahedral  Elements 

A  typical  tetrahedral  element  is  shown  in  Figure  16.  It  is  geometrically  defined  by  nodes  i,  j,  m, 
and  p.  The  formulation  is  based  on  nodes  i,  j,  and  m  being  positioned  in  a  counterclockwise 
manner  when  viewed  from  node  p.  The  coordinates  of  node  i  are  designated  Xi,  yi,  zj,  and  the 
corresponding  velocities  are  designated  Uj,Vi,Wj  (which  are  identical  to  Xj,yi,Zj  used 

previously).  This  algorithm  is  reported  in  References  2, 4, 14,  and  23. 
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The  algorithm  that  follows  is  for  a  single  tetrahedral  element.  The  arrangement  of  the  elements, 
however,  can  have  a  significant  effect  on  the  computed  response  (References  16  and  18).  There 
are  also  some  mixed  (average  volumetric  strain,  or  average  pressure)  algorithms  that  provide 
improved  acciuucy  for  various  arrangements  of  elements  (Reference  16). 

The  mass  at  node  i,  for  an  individual  tetrahedral  element,  is  simply 

Mi=p„V„/4  (1) 


where  po  is  the  initial  density  of  the  material  and  Vo  is  the  initial  volume  of  the  element.  The 
general  expression  for  the  volume  is 


V  =  - 


1  Xj  y;  Zi 

1  X.  y.  Zj 

1  y„  z„ 
1  X„  v„  z„ 


(2) 


When  the  elements  are  incorporated  into  an  assemblage  of  elements,  then  the  total  mass  at  node  i 
contains  individual  element  masses  from  all  elements  that  contain  that  node. 

(3) 


The  volumetric  strain  and  strain  rate  are  obtained  in  the  same  manner  as  used  for  the  ID  and  2D 
elements. 

Ev=V/V„-l  (4) 

=  (e^''"-et)/At  (5) 

In  Equation  4,  V  and  Vo  are  the  current  and  initial  element  volumes;  and  in  Equation  5,  and 
are  the  volumetric  strains  at  time  t  +  At  and  time  =  t.  The  volumetric  strains  and  strain  rates 
are  based  on  the  initial  configuration,  as  opposed  to  the  shear  and  deviator  strain  rates  which  are 
based  on  the  current  configuration. 
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The  strain  rates  are  obtained  from  the  current  geometry  of  the  element  and  the  velocities  of  the 
nodes.  If  it  is  assumed  the  velocities  vary  linearly  between  the  nodes,  the  x,  y,  and  z  velocities 


)  within  each  element  can  be  expressed  as 

u  =  a,  +a2X+a3y+a4Z 

(6) 

V  =  ttj  -I-  agX + a^y + ttgZ 

(7) 

w  =  a,  +  ttjpX  ■+■  a„y + a,2z 

(8) 

where  ai  . . .  an  are  geometry  and  velocity-dependent  constants  for  each  element.  It  is  possible 
to  solve  for  ai  ...  04  by  substituting  the  x  velocities  and  coordinates  of  the  four  nodes  into 
Equation  6.  This  gives  four  equations  and  four  unknowns  so  that  the  constants  (ai  ...  04)  can  be 
evaluated. 


Equation  6  can  then  be  expressed  in  terms  of  element  geometry  and  nodal  velocities. 

u  =  ^[(ai  +biX+Ciy-i-diz)ui  +(aj  -HbjX+c.y+djz)uj 

+  (am  +  b„x + c„y + d„z)u„  +  (ap  +  bpX + Cpy + dpz)up  ] 


where  V  is  the  volume  from  Equation  2  and  the  geometry  constants  have  the  form 


ai  = 


yj 

Ym 


Zj 


Xp  Yp  Zp 


(10) 


bi=- 


1  Yi 

1  Ym  Zm 

1  Yp  Zp 


(11) 


Ci  = 


1  X; 


Zj 


(12) 
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(13) 


1  Xj  Yj 

1  X„  Vm 

1  Xp  Yp 


The  remaining  geomehy-dependent  constants  (aj,  bj,  cj,  dj,  etc.)  are  obtained  bY  a  sYstematic 
interchange  of  signs  and  subscripts.  The  y  and  z  velocities  in  Equations  7  and  8  are  obtained  in  a 
similar  manner  and  are  identical  to  Equation  6  except  the  x  velocities  at  the  four  nodes  are 
replaced  bY  the  y  and  z  velocities. 


Another  geomehy  constant  is  the  altitude  of  the  tetrahedral  element.  The  altitude  from  node  i  to 
the  plane  defined  bY  the  other  three  nodes  of  the  tetrahedron  is 


h;  = 


6V 


Vbf  +cf  +df 


(14) 


The  other  altitudes  (hj,  hm,  hp)  are  obtained  bY  appropriate^  changing  the  subscripts  of  the 
geometty-dependent  constants.  The  minimum  of  the  four  altitudes  is  designated  as  h,  and  it  is 
used  for  the  characteristic  length  in  Equation  10  of  subsection  3.8.1,  for  the  integration  time 
increment. 

After  the  velocities  are  obtained,  it  is  possible  to  determine  the  normal  strain  rates  (e^,  Ey,  e^) , 
the  shear  strain  rates  (y^y,  y„,  Yy^)  and  the  spin  rates  (©^,0)y,0)J  of  the  element. 


(15) 

(16) 

(17) 

9u  dv 
dy  9x 

(18) 
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3u  3w 


(19) 


9v  3w 


(20) 


1 

CO.  =  2 


( dw 

9y  dz 


(21) 


ifdii  Bw^ 


(22) 


“■=i 


dv  9u 
dx  By 


(23) 


It  can  be  seen  that  Equations  15-23  are  derivatives  of  linear  functions  and  therefore  give  constant 
values  within  the  element. 


The  equivalent  strain  rate  (e)  and  deviator  strain  rates 


e.,e. 


are  determined  in  an  identical 


manner  as  shown  previously  for  ID  and  2D  geometry. 


(24) 


The  constitutive  models  also  require  deviator  strain  rates,  which  are  expressed  as 


®x  ^ave 


®y  ^ave 


(25) 

(26) 


®z  ^ave 


(27) 


where  =(e,j+ey+e^)/3.  Note  that  the  sum  of  the  deviator  strain  rates  is  e^  +  By  +  e^  =  0 . 


B  t19616.doc 


81 


The  stresses  in  the  elements  are  determined  from  the  strains,  strain  rates,  temperatures,  pressures, 
internal  energies,  and  material  constants.  This  is  identical  to  the  ID  and  2D  algorithms  for  solid 
elements,  except  that  all  three  shear  stresses  are  present.  The  three  normal  stresses  are  generally 
expressed  as 


o.=s,-(P  +  Q) 

(28) 

o,  =s, -(p  +  Q) 

(29) 

<j.  =  s,-(p  +  Q) 

(30) 

where  Sx,  Xy,  and  Sz  are  the  normal  deviator  stresses,  P  is  the  hydrostatic  pressure  and  Q  is  the 
artificial  viscosity.  These  are  described  in  detail  in  Section  6. 


Trial  values  of  the  deviator  stresses  and  shear  stresses  at  time  =  t  +  At  are 

sr=sL+2Ge,At  +  As, 

s;^^=s;+2GeyAt  +  ASy 

s*-"^' =s‘,+2Ge,At  +  As, 

<r=<y  +  GYxyAt  +  Ax,y 
='cL+GYxzAt  +  AT,, 

'^T  ='»^yz  +  GYyzAt  +  AXy, 


(31) 

(32) 

(33) 

(34) 

(35) 

(36) 


In  Equation  3 1  the  first  term  (s* )  is  the  normal  stress  at  the  previous  time  and  the  second  term 
(2Ge^At)  is  the  incremental  stress  due  to  the  incremental  strain  (e,,At)  during  that  time 

increment,  where  G  is  the  elastic  shear  modulus  and  At  is  the  integration  time  increment.  The 
third  term  (as,j)  is  due  to  shear  stresses  from  the  previous  time  increment,  which  now  act  as 

normal  stresses  due  to  the  new  orientation  of  the  element  caused  by  an  incremental  rotation 
^(OyAt,  (OzAtj  during  the  time  increment.  The  remaining  normal  stresses  and  shear  stresses  have 

a  similar  form. 


B  t19616.doc 


82 


The  correction  terms  for  element  rotations  are 


As,  =2(©y<-®,<jAt 

(37) 

ASy=2(Q)Xy-®x'c;z)At 

(38) 

As,  =2  ((0,1;,-©/,,)  At 

(39) 

=[©,  (s* -s^j  +  mXz-^x'cL 

|At 

(40) 

Ax„  =  [©y  (s;  -  sO  +  co,t‘ y  -  ©,t;. 

At 

(41) 

AV  =[a),  (s;  -s‘)  +  ©,t;^  -tOy'Cxy 

At 

(42) 

Equations  3 1-36  assume  an  elastic  response  of  the  material.  If  the  strength  of  the  material  is 
exceeded,  then  plastic  flow  (or  fracture)  will  occur.  The  Von  Mises  yield  criterion  is  used  to 
determine  an  equivalent  stress,  o ,  that  can  be  compared  to  the  uniaxial  tensile  (or  compressive) 
strength  of  the  material.  The  general  form  of  the  equivalent  stress  is 

^  =  +(<yx  -  +  K  -  <7zr  +Ty,')]  (43) 

Using  deviator  stresses  (instead  of  total  stresses).  Equation  43  can  be  rewritten  as 

a  =  ^|(s^  +  sj  +  ) + 3(x',y  +  tL  +  xl,)  (44) 

If  a  is  not  greater  than  the  equivalent  tensile  strength  of  the  material,  ct,  the  final  deviator  and 
shear  stresses  are  as  given  in  Equations  31—36.  If  o  is  greater  than  ct,  then  the  stresses  in 
Equations  3 1-36  are  multiplied  by  the  factor  (a  /  a) .  When  the  reduced  deviator  and  shear 
stresses  are  put  into  Equation  44,  the  result  is  always  a  =  a .  This  is  known  as  the  radial  return 
algorithm.  The  various  material  strength  models  for  ct  are  presented  in  Section  6. 


B_t19616.doc 


83 


During  plastic  flow,  it  is  sometimes  necessary  to  determine  the  equivalent  plastic  strain  for  strain 
hardening  effects  on  the  strength  of  the  material  or  to  determine  if  the  material  has  failed.  The 
first  step  in  this  process  is  to  adjust  the  total  strain  rates  to  plastic  strain  rates  by  subtracting  out 
the  elastic  portion  of  the  strain  rates. 


(45) 

4;=4,-(s"*-s;-As,)/2GAt 

(46) 

4J=4.-(s“’-s;-AsJ/2GAt 

(47) 

(48) 

(49) 

Y^=Y„-(^r-';x-At„)/GAt 

(50) 

The  general  expression  for  the  equivalent  plastic  strain  rate  is 


V9 


(4;  -4;)'  +(4;  -4;)’ +(4;  -4;)+-(y:/  +rJ  +Y^") 


(51) 


The  equivalent  plastic  strain,  £p ,  is  then  obtained  by  integrating  £p  with  respect  to  time. 


£;""‘=£;+£pAt 


(52) 


After  the  element  stresses  are  determined,  the  concentrated  nodal  forces  can  be  obtained.  These 
forces  are  statically  equivalent  to  the  distributed  stresses  within  the  element.  They  are  dependent 
on  the  displaced  element  geometry  and  the  magnitude  of  the  stresses.  The  forces  in  the  x,  y  and  z 
directions  at  node  i  of  an  element  are  given  by 


+  C:T 


xy 


(53) 
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(54) 


F;=-^(ci<?y+bi'Cxy+diV) 

Fi  =  -^(dia,+CiXy,+biT„)  (55) 


The  geometry-dependent  constants  (bj,  Cj,  di)  are  again  identical  to  those  used  for  calculation  of 
the  strain  rates  and  altitudes.  The  forces  at  the  other  nodes  are  readily  obtained  by  changing 
subscripts. 

Note  that  -t-  +  F”  +  F^**  s  0  for  each  element,  and  this  ensures  equilibrium  in  the  system.  A 

similar  equilibrium  exists  in  the  y  and  z  directions.  This  results  in  conservation  of  momentum  if 
there  are  no  external  forces  or  restraints.  The  final  net  forces  on  node  i  are 


F'=SFi 


(56) 


(57) 


Fj=SFi 


(58) 


It  is  also  possible  to  add  external  forces  through  applied  pressures. 


The  heat  conduction  algorithm  can  be  obtained  straightforwardly  fi'om  the  ID  heat  conduction 
algorithm  (subsection  3.1.3)  and  the  strain  rate  equations  in  this  subsection.  The  2D  algorithm  is 
reported  in  Reference  7.  Substituting  temperatures  for  velocities  in  Equation  9  gives 


^  +  biX+ c^y + diz)T.  +  (aj  +  bjX+ c^y + djz)Tj 

+  (a™  +  b„.x + c,y  +  d„z)T„  +  (a^  -h  b^x + Cpy + dpz)Tp] 


(59) 


Other  than  the  nodal  temperatures  (Tj,  Tj,  Tm,  Tp)  all  of  the  other  terms  are  as  described  for 
Equation  9. 

The  instantaneous  heat  flows  in  the  x,  y,  z  directions  can  then  be  obtained  from 
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q.  =  -k|^  =  -k(b,T,  +b^Tj  +b.T.  +b,Tj/6V 

(60) 

qy  =  =  -k(ciTi  +CjTj  +  c„T„  +  CpTp)  /  6V 

(61) 

q  =-k— =  -kfdT +d  T +d  T  +dTl/6V 

Mz  ^  1  j  j  m  m  p  p  / 

(62) 

The  incremental  increase  in  thermal  energy  at  the  nodes  can  be  obtained  by  integrating  the  heat 

flow  with  respect  to  area  and  time. 

AQi  =|(qxbi  +qyCi  +q,di)At  +  AQ,  /4 

(63) 

Note  that  the  internal  energy  generated  in  the  element  during  the  previous  time  increment,  AQe,  is 
distributed  equally  (AQe/4)  to  all  four  nodes. 

The  remaining  equations  are  similar  to  those  in  ID  and  2D.  The  updated  temperatures  of  the 
nodes  have  the  form 

^Q.  / (64) 

where  T/*^*  and  T/  are  the  temperatures  of  the  nodes  at  times  t+At  and  t,  ZAQj  is  the  sum  of  the 
incremental  heat  contributed  by  all  elements  that  contain  node  i,  Mj  is  the  total  mass  of  node  i, 

and  Cpi  is  the  specific  heat  of  node  i. 

The  internal  energy  in  an  element  can  be  used  to  compute  element  pressures.  To  account  for  the 
flow  of  internal  energy  through  the  grid,  the  element  temperature  is  assumed  to  be  the  average  of 
the  nodal  temperatures. 

T  =  (Ti+Tj  +  T,+Tp)/4  (65) 

The  internal  energy  (per  initial  volume)  is  given  by 

E.=(T-T,)p.Op  (66) 


B  t19616.doc 


86 


where  To  is  the  initial  temperature  (where  Eg  =  0),  po  is  the  initial  density,  and  Cp  is  the  specific 
heat  of  the  element  material. 

The  integration  time  increment  must  also  be  bounded  to  ensure  that  the  computations  remain 
stable  for  heat  conduction  (References  6  and  7).  The  heat  conduction  portion  requires 

At<pcXin/4k  (67) 

where  hmin  is  the  minimum  altitude  of  the  element  (as  provided  in  Equation  14)  and  the  other 
terms  recently  have  been  defined.  This  is  analogous  to  the  time  increment  restriction  for  wave 
propagation  in  Equation  10  of  subsection  3.8.1.  Unless  hmin  is  very  small,  the  wave  propagation 
restriction  is  much  more  severe  than  the  heat  conduction  restriction. 

3.8.5  Brick  Elements 

The  3D  brick  element  is  shown  in  Figiure  19.  This  element  was  developed  by  Flanagan  and 
Belytschko  (Reference  26)  and  the  nodal  numbering  in  Figure  19  is  as  used  in  their  paper.  The 
EPIC  code  uses  a  different  nodal  numbering  arrangement  for  the  user,  and  the  transformation  is 
made  within  the  code.  In  Figure  19,  nodes  5, 6, 7,  and  8  are  clockwise  when  viewed  from  node 
1.  The  coordinates  of  node  1  are  xi,  yi,  zi,  and  the  corresponding  velocities  are  u, ,  v, ,  w, 
(which  are  identical  to  x, ,  y , ,  z, ). 

The  mass  at  node  i,  for  an  individual  brick  element,  is  simply 

Mi=p„V„/8  (1) 

where  p,,  is  the  initial  density  of  the  material  and  Vo  is  the  initial  volume  of  the  element.  The 
volume  is  a  complex  expression  that  is  presented  later. 

The  strain  rates  are  obtained  fi'om  the  current  geometry  of  the  element  and  the  velocities  of  the 
nodes.  Note  that  the  geometry  is  expressed  in  a  reference  domain  coordinate  system  as  a  unit 
cube,  as  shown  in  the  upper  right  portion  of  Figure  19.  The  velocity  in  the  x  direction  is 

U  =  (t)iU,  +(^>2^2  +<t>3U3  +<>4^4  +<(>5^5  +<1)6^6  +<l)7il7  (2) 
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Reference  Domain 


Physical  Domain 


X 

28.  vsd 
T19616 


Hourglass  Modes 

Z  z 


Mode  3  (X)  Mode  4  (X) 

29.vsd 

T19616 


Figure  19.  Description  of  the  3D  Brick  Element 
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where  the  shape  functions  are 


(3a) 

(3b) 

(3c) 

(3d) 

(3e) 

(3f) 

(3g) 

(3h) 

where  and  r\  are  the  coordinates  of  the  reference  domain  as  shown  in  Figure  19.  The  y  and  z 
velocities  (v  and  w)  have  the  same  form  as  Equation  2.  The  x  nodal  velocities  (u,  ...  lig)  are 

simply  replaced  by  the  y  nodal  velocities  (vj ...  Vg)  and  the  z  nodal  velocities  (wj ...  Wg). 

For  the  physical  domain  it  is  necessary  to  determine  the  [B]  matrix  for  the  geometry  constants. 


B.2 

Bx3 

B.4 

B.5 

B,5 

B.7 

B.S 

[B]  = 

B,2 

By3 

By4 

ByS 

By6 

By7 

ByS 

B.1 

B.2 

B.3 

B.4 

B.5 

B., 

B.7 

Bz3. 

B  t19616.doc 


89 


where 


B„=y2[(z6-Z3)  +  (z5-Z4)]+y3(z2-Z4)  +  y4[(z3-Z8)  +  (z2-Z5)] 

+  y5[(z8-Z6)+(z4  -Z2)]+y6(z5-Z2)+y8(z4  -Z5) 

B,2  =  ysRzv  -  Z4  )  +  (z6  -  Z, )]+  y4(z3  -  Z, )+  y,[(z4  -Z,)  +  (z,-  zj] 
+  ye  [(zj  -  Z7  )  +  (z,  -  Z3 )]  +  y-,  (z,  -  Z3 )  +  Ys  (z,  -  z,  ) 

B,3  =  y4 [(zg  - z, )+ (z^  - Z2 )]+  y, (z,  -z^)+y, [(z,  -z,)+(z,- z, )] 
+y7[(z6  - Z8)+ (z2  -  Z4)]+y8(z7  -  z,)+  y^iz^  - z,) 

Bx4  =  yi[(z5  -Z2)  +  (z8  -Z3)]+y2(zi  -Z3)+y3[(z2  -Z7)  +  (z,  -zj] 

+  y*  [(z^  -  Z5 )  +  (z3  -  z, )] + Ys  (zg  -  z, ) + y 7  (z3  -  Zg ) 

Bx5  =  yg  [(z4  -  Z5 ) + (z,  -  Zg )] + y7  (zg  -  Zg ) + y g  [(z^  -  z^ ) + (zg  -  z, )] 
+  y I  [(z2  -  Z4  )  +  (Zfi  -  Zg  )] + y4  (z,  -  Zg  )  +  y2  (zg  -  z, ) 

Bx6  =  ys^Zi  -Zg)+(z2  -Z2)]+yg(z5  -Z7)+y2[(zg-Z3)+(z5  -Z2)] 

■*"y2[(z3  ~Z,)+(z7  “Z5)j  +  yj(z2  — Z5)+y3(z7  — Z2) 

Bx7  =  y6[(z2  -Z5)+(z3  -Zg)]+y5(zg  -Zg)+yg[(z5  -Z4)+(zg  -Z3)] 
+  y3[(z4  -  Z2)+  (Z8  -  Zg)]+ y2  (z3  -  Zg)+  y4(zg  -  Z3) 

Bxg  =  y7[(z3  -  Z6)  +  (Z4  -  Z5)]+  yg(z7  -  Z5)  +  y5[(zg  -Z,)  +  (z2  -  Z4)] 

+  y^  [(z,  -  Z3 )  +  (z5  -  Zy )] + y  3  (z4  -  z, ) + y ,  (z^  -  Z4 ) 


(5a) 

(5b) 

(5c) 

(5d) 

(5e) 

(5f) 

(5g) 

(5h) 


For  By, . .  .Byg  replace  y  and  z  with  z  and  x,  and  for  B^, . . .  B^^g  replace  y  and  z  with  x  and  y. 


The  volume  can  now  be  determined  using  any  of  the  following  three  expressions. 
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(6a) 


V  =  +  6,2X2  +  6,3X3  +  6,4X4  +  6,5X5  +  6,5X,  +  6,2X2  +  6,8X  J 

V  =  ]^{By,y,  +  By2y2  +  By3y3  +  6y4y4  +  6y5y5  +  +  6y8y8}  (6b) 

V  =  ^{b„z,  +  6,2Z2  +  6,323  +  6,424  +  6,525  +  6,525  +  6,222  +  6,328}  (6c) 


The  final  geometric  term  is  the  characteristic  element  length,  h,  that  is  used  for  the  artificial 
viscosity  and  the  integration  time  increment. 


h  = 


(7) 


where  V  is  the  element  volume  and 

6,^  =  6„^  +  6,2^  +  6,3^  +  6,4^  +  6,5^  +  6,5^  +  6,2^  +  6,3^ 

+  6,,^  +  6,2'  +  B,3^  +  B,4^  +  6^5^  +  6^5^  +  6,2^  +  6,3^  (8) 

+  6„^  +  6,2^  +  6,3^  +  6,4^  +  6,5^  +  6,5^  +  6,2^  +  6,3^ 


The  three  normal  strain  rates  (e,  ,  £y ,  e, )  are  obtained  by  taking  derivatives  of  the  velocity 
fields  in  the  physical  domain. 


Ex  =  +  ®x2U2  +  B,3U3  +  6,4114  +  6,5*5  +  B,5U5  +  6,2*2  +  ]  (9) 


dv  1 


ay  12V 


[6y,V,  +  6y2V2  +  6^3*3  +  6^4*4  +  6^5*5  +  6^5*5  +  6^2  V2  +  6^3*3]  (10) 


e.=~  =  t:^[B„  w,  +  6,2  W2  +  6,3W3  +  6,4  W4  +  6,5  W5  +  6,5W5 


92  12V 


+  6,2  W2 +6,3^3] 


(11) 
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The  X  velocity  in  Equation  9  is  expressed  as  a  function  of  the  reference  domain  coordinates  in 

f  96,  9(|)o 

Equation  2.  The  derivatives  of  the  terms  in  Equation  2  are  replaced  by  the 

physical  domain  geometry  constants  in  Equation  9  •  For  the  y  direction. 


9<t)i  9(l)g 


B 


B 


9(|),  9(1)8 


9y  9y 
replaced  by 


are  replaced  by  — — . . .  — —  in  Equation  1 0  and  for  the  z  direction,  -:r~ . . .  ~r—  are 
12V  12V  9z  9z 


B 


zl 


B 


z8 


12V"  12V 


in  Equation  1 1. 


The  three  shear  strain  rates  y^,  Yj^)  and  the  three  rotational  spin  rates  (d)^,  (by,  (b^)  are 
determined  in  a  similar  manner. 


9u  9v 


(12) 


9u  9w 


(13) 


9v  9w 


(14) 


CO  =  — 
*  2 


1  f  9w  ^ 
9y  9z^ 


(15) 


_  1 1^911  9w' 
2I9Z  9x, 


(16) 


ir9v  9u^ 


2I9X  9y 


(17) 


The  stresses  are  determined  with  exactly  the  same  procedures  used  for  the  tetrahedral  elements  in 
subsection  3.8.4. 


The  X  direction  forces  on  each  of  the  eight  nodes  are  as  follows: 
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1 

F  =  - — 
12 

Bx,Cfx+By,T,y+B„x,„]-F,f 

(18) 

1 

“  “  12 

®x2^x  ■^®y2'^xy  ®z2'^xz  ]  ~  Px2 

(19) 

1 

““l2 

Bxs^^x  +  Bys'Cxy  +  B^jX,^]  -  F,f 

(20) 

1 

F  _ _ 

12 

[Bx4*^x  By4X,(y  +  B^4X,(2]  — F,(4 

(21) 

1 

F  =  - — 
12 

BxsCJx+Bys'Cxy+B^x^J-F.f 

(22) 

„  1 

12 

Bx6<^X  +  By6X,jy  +  B^gX,^]  -  F^’f 

(23) 

^  1 
12 

[Bx7^X  "^ByyX^jy  +B^7X,j2]“Bx7 

(24) 

1 

”  12 

Bx8*^x  Byg'txy  BzgX,^^  j  —  F^g 

(25) 

where  B^, . . .  B^g  are  the  geometry  constants,  CTx  is  the  normal  stress  in  the  x  direction  and  x 

xy 

and  are  shear  stresses.  The  hourglass  forces  (F”°...Fjg°)  are  required  to  resist  hourglassing 
deformations  and  they  are  presented  later.  The  y-direction  forces  (pyl-’-fys)  are  obtained  by 
replacing  a^,  andx^  by  x^y,  Oy,  andXy^  respectively,  and  the  z-direction  forces 
(F^p.-F^g)  are  obtained  by  replacing  x^y,  andx^  with  x^,  Xy^,  and  a  respectively. 

A  3D  brick  element  can  experience  hourglass  deformations  in  a  manner  similar  to  a  2D  quad 
element.  In  3D,  however,  it  is  much  more  complex.  The  lower  portion  of  Figure  19  shows  four 
homglass  modes  in  the  x  direction.  The  same  four  modes  occur  in  the  y  direction  and  in  the  z 
direction,  for  a  total  of  12  hourglass  modes. 

The  following  four  steps  are  required  for  the  hourglass  algorithm: 
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.  Compute  the  hourglass  shape  vector  (y) 

•  Compute  the  hourglass  modal  velocity  (“i) 

•  Compute  the  hourglass  resistance  stress  (Q) 

•  Compute  the  hourglass  resistance  force 


The  hourglass  shape  vector  for  mode  1  is 

Tn  =  +  (26a) 

Yi!=l-;^[B.!>‘+B,2y+B.iz]  (26b) 

Y„  =  -1-^[B.35E+B„y  +  B^z]  (26c) 

Yi4=-1— [|^[B,<x+B„y  +  B,,z]  (26d) 

Yi.=-1-^[B.6X+B„y+B,.z]  (26f) 

Y.>=l-f^[B,.x+B„y  +  B.,z]  (26g) 

Yit  =  l-J^[B..x+B„y+Baz]  (26h) 

where  . . .  B^g  are  the  geometry  constants  and 

X  =  X,+X2-X3-X4-X5-X6+X7+X8  (27) 

y  =  yi+y2-y3-y4-y5-y6+y7+y8 
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Z  =  Zj  +  Z2  *“  Z3  “  +  Z-7  +  Zg 


(29) 


The  shape  vectors  (y,, . .  .y,g )  in  Equation  26  are  not  related  to  the  shear  strain  rates 

(txyj  Yxz>  Yyz)  ^  Equations  10-12.  The  shape  vectors  for  modes  2, 3,  and  4  (721. ..728,  Y3i...y38, 


and  741 . .  .748)  are  obtained  in  a  similar  manner. 

The  hourglass  modal  velocities  for  the  fom  x  modes  in  Figure  19  are 

4x1  ~  Yi2'^2  Yi3^3  Yi4^4  Yi5^5  ■*"  Yi6^6  ■*"  Yl?'^?  Ylg^s] 

4x2  =;^[Y2lU,+Y22U2+Y23U3+Y24U4+Y25U5+Y26i»6+Y27U7+Y28Ug]  (31) 

4x3  ~  [Y31^1  Y32^2  Y33'^3  Y34^4  Y35'^5  Y36^6  Y37^7  "*■  Y38^8]  (32) 

4x4  =^[Y4iU.  +Y42U2+Y43U3+Y44U4+Y45U5+Y46U6+Y47U7+Y48U8]  (33) 


The  y  modal  velocities  (qy, . . .  qy4 )  are  obtained  by  replacing  u  with  v ,  and  the  z  modal 
velocities  (qzi...qz4)  are  obtained  by  replacing  u  with  w. 

The  resisting  stress  for  the  first  x  mode  is 

Qxi  ~  ^hP^s4xiV®x  (34) 

where  Ch  is  an  input  hourglass  viscosity  coefficient,  Cs  is  the  sound  velocity  of  the  material,  p  is 
the  density  of  the  material,  and 

Bx  =  [BxI  +  +  B^3  +  B^4  +  B^5  +  B^6  +  B^7  +  B^g  ]  (35) 

The  resistance  stresses  for  the  other  three  x  modes  (q^j  >  Qx3  »  Qx4 )  obtained  by  replacing 
q,;,  by  the  other  modal  velocities  (q^j »  4x3  >  4*4 )  •  The  y  and  z  modes  are  obtained  in  a  similar 
manner. 
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Now  the  forces  to  resist  hourglassing  can  be  determined  and  substituted  into  Equations  18-25. 


“■^^[QxiYiI  ■*'Qx2Y21  ■*'Qx3Y31  ■*‘Qx4Y4i] 

(36) 

~  ^[QxiYi2  ■*'Qx2Y22  ■*'Qx3Y32  ■^Qx4Y42] 

(37) 

~  ^  Qx2Y23  Qx3Y33  ■*‘Qx4Y43] 

(38) 

“■^^[QxiYm  ■^Qx2Y24  ■*"Qx3Y34  +QX4Y44] 

(39) 

=  ■*'Qx2Y25  ■•■Qx3Y35  ■*'Qx4Y45] 

(40) 

—  ^^[QxiYi6  ■*‘Qx2Y26  ■*‘Qx3Y36  ■*’Qx4Y46] 

(41) 

~  ^^[QxiYi7  Qx2Y 27  Qx3Y 37  Qx4  Y 47  ] 

(42) 

=  ^[QxiYi8  ■*‘Qx2Y28  ■*"Qx3Y38  ■*‘Qx4Y48] 

(43) 

For  Ff  ...Ff  replace  Q„...Q,4  with  Qy,...Qy4,  and  for  Ff...Ff  replace  Q„...Q,4  with 
Qzr--Qz4  • 


Both  the  homglass  forces  and  the  internal  stress  forces  are  then  added  together  as  shown  in 
Equations  18-25. 


The  heat  conduction  algorithm  follows  the  same  pattern  as  used  for  the  other  elements.  The  heat 
flow  in  the  three  directions  is 
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(44) 


®x3"^3  ®x4"r4  ®x5"^5  ■*"  ®x6"^6  ■*■  ®x7^7  "*■  ®x8 


Tg]/12V 


qy  =-k 


8y 


•  “k[By,T,  +  By^Tj  +  ByjTj  +  B  y  4  T4  +  B  y  J  +  B  y  ^  +  B  y  ^  +  B  y  g  ig  ]  /  ^  V 


(45) 


-f 


~  ■*■  ®22^2  ®z3"^3  ■*■  ®Z4"^4  ■*"  ^zS^^S  ■*■  ^26^6  +  Bz7"^7  "*■  ^zS^^s]  ^ 


(46) 


The  increase  in  thermal  energy  at  node  1  is  AQ,  =  ^(B^,q^  +  By,qy  +  B^,q^)At  +  AQ^  /  8 ,  where 

Bxi,  Byi,  and  Bzi  are  as  described  previously,  At  is  the  integration  time  increment,  and  AQe  is  the 
internal  energy  generated  by  the  element  during  the  previous  time  increment.  The  temperature 
and  internal  energy  update  are  similar  to  those  described  for  the  tetrahedral  elements  in 
subsection  3.8.4. 


3.8.6  Bar  Elements 

The  3D  bar  element  is  a  simple  element  that  is  incompressible  and  carries  only  axial  loads.  It  is 
described  in  Reference  27.  A  summary  of  the  3D  bar  element  algorithm  is  shown  in  Figure  20. 
The  bar  is  geometrically  determined  by  nodes  i  and  j  (whose  positions  define  the  length,  i ),  and 
the  cross-sectional  area,  A.  By  assiiming  the  material  is  incompressible,  the  area  is  determined 
firom 


(1) 
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where  V,,,  A„,  and  are  the  initial  volume,  initial  cross-sectional  area,  and  initial  length, 
respectively.  The  mass  is  equally  distributed  to  the  two  nodes. 


Mi=M-p„V,/2  (2) 

where  po  and  Vo  are  the  initial  density  and  volume. 

If  the  z'  axis  is  assumed  to  coincide  with  the  axis  of  the  bar,  the  axial  strain  rate  is  given  by 

(3) 

where  Vi  and  Vj  are  the  velocity  components  of  nodes  i  and  j,  along  the  axis  of  the  rod. 
Specifically 

Vj  =(uif„ -t-Vj^y -l-w/ J  (4) 

where  Uj ,  Vj ,  and  Wj  are  the  x,  y,  and  z  components  of  the  velocity  at  node  i  and  £^,  £y,  £^  are 
the  direction  cosines  of  the  bar. 

Incompressibility  requires  that 

e,.+£y.-he,,=0  (5) 

where  and  Ey-  are  the  strain  rates  normal  to  the  longitudinal  axis  of  the  bar.  This  gives 

Ex' 2  (6) 

with  the  strain  rates  defined,  the  deviator  stresses  can  be  determined  by  the  standard  method. 

The  net  stresses  are  a  function  of  both  the  deviator  stresses  fs^.,s„.,s^J  and  the  pressure  (P). 

Or  =  s..-P  (7) 

d..=s..-P  =  0  (8) 
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Because  the  net  normal  stresses  must  vanish  =  o) ,  the  pressure  is  defined  as 

P  =  s^-=Sy.  (10) 

The  net  nodal  forces  act  along  the  axis  of  the  bar  and  are  determined  from 

Fi=Fj=a,-A  (11) 


The  components  in  the  x,  y,  and  z  directions  are  given  by 


F' 

^  X 


X 


(12) 


F'  =  F  ^ 

y  1  y 

F‘  =  F^ 

where  are  the  direction  cosines  of  the  bar. 


(13) 

(14) 


The  heat  conduction  algorithm  is  identical  to  the  ID  heat  conduction  algorithm  in  subsection 
3.1.3,  except  the  area  for  the  3D  bar  is  given  by  Equation  1 . 

3.8.7  Bending  Shell  Elements 

A  3D  bending  shell  element  is  shown  in  Figures  21-25.  The  concepts  for  the  algorithm  are 
presented  in  Reference  20.  As  was  the  case  for  the  2D  bending  shell  element,  described  in 
subsection  3.4.6,  the  3D  bending  shell  is  represented  by  three  or  five  layered  elements  that  share 
the  same  nodes.  Each  of  the  three  or  five  layered  elements  is  a  triangle  defined  by  nodes  i,  j,  k. 
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z 


44.VSCI 
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Figure  23.  Description  of  the  Mass  Distribution  for  the  3D  Bending  Shell  Element 
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z 


Sign  Convention 

Positive  stress  on  positive  face  positive  sign 
Positive  stress  on  negative  face  negative  sign 
Positive  moment  about  axis  ->  right  hand  rule 


z 


45.vsd  CT2 

T19616 


Figure  24.  Description  of  the  Stresses  for  the  3D  Bending  Shell  Element 
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46.vsd 

T19616 


Figure  25.  Description  of  the  Forces  and  Moments  for  the  3D  Bending  Shell  Element 

The  first  step  is  to  perform  a  coordinate  transformation  into  a  local  coordinate  system,  as  shown 
in  Figure  21 .  Each  triangle  is  put  into  the  x'  -  z'  plane  at  y'  =  0 ,  where  node  i  is  at 
x'  =  y'  =  z'  =  0 ,  and  node  j  is  on  the  positive  x'  axis  at  y'  =  z'  =  0 .  The  transformations  are  not 
presented  herein.  They  are  partially  provided  in  Reference  27. 


Unless  noted  otherwise,  the  remainder  of  this  subsection  is  presented  in  the  local  coordinate 
system,  and  the  x',  y',  z'  coordinates  will  be  represented  as  x,  y,  z  for  simplicity. 
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Bending  moments  in  the  shell  element  are  induced  by  allowing  the  in-plane  stresses  in  the  outer 
layers  to  be  different  from  one  another.  This  is  accomplished  by  allowing  the  in-plane  strain 
rates  to  be  different  from  one  another.  These  strain  rates  are  dependent  on  the  nodal  velocities, 
as  well  as  the  nodal  rotational  rates  and  offset  from  the  center  layer,  as  shown  in  Figure  22. 

The  masses  at  the  three  nodes  of  an  element  are  equal. 

Mi=Mj  =  Mk=p„V„/3  (1) 


where  po  is  the  initial  density  and  Vo  is  the  initial  volume.  The  initial  volume  is  simply  Vo  - 
AoTo,  where  Ao  is  the  initial  area  and  To  is  the  initial  total  thickness  (of  all  layers). 

The  mass  moments  of  inertia  are  required  for  the  rotational  rates  at  the  nodes.  Figure  23  shows 
how  the  mass  is  distributed  to  the  three  nodes  for  the  inertia  computations.  The  mass  moments 
of  inertia  for  node  i  are  as  follows: 


4=pbhT 


23h^  T^l 
2592'*’ 72  J 


(2) 


lL=lL=-pbhT 


14bh-i-23hxK 


where  p  is  the  current  density,  T  is  the  current  total  thickness  of  the  element  (all  layers),  and  b,  h, 
and  Xk  are  geometric  distances  provided  in  Figure  23. 
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The  moments  of  inertia  for  node  j  are 


i  f23h^  tM 

Ixx  -  Pbhlj  +  72  I 


(8) 


4=pbhP 


222b"  +69h^  -222bx,  +69x,^ 
7776 


(9) 


lL=pbhT 


[(222b"-222bXk+69Xk^)  t" 
7776 


(10) 


lL=lL  =  pbhT^ 


37bh-23hx;, 

2592 


(11) 


liy=lix=0  (12) 

liz=liy=0  (13) 


and  the  moments  of  inertia  for  node  k  are 


37h"  T" 


(14) 


4=pbhT 


46b"  +  148h"  -  148bx^  +  148x^ 
5184 


(15) 


lL=pbhT 


[(46b"-148bXk+148x;,")  T" 


5184 


72 


(16) 


l!L=lL=pbhT 


37h(b-2xJ] 


2592 


(17) 


1“  =1*'  =0 

xy  yx 


(18) 
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(19) 


l'‘  =l‘‘  =0 

Note  that  the  current  (not  initial)  geometry  must  be  used  and  this  requires  Ixx  . .  .Izz  to  be  updated 
every  cycle. 

In  Figure  22  the  three  comer  velocities  of  layer  1  are  identical  to  the  nodal  velocities,  because  the 
nodes  are  positioned  in  the  plane  of  layer  1.  For  nomenclature  consistency,  the  nodal  x  velocities 
at  node  i  are  designated  Uj ,  Ujj ,  Uj3 ,  the  y  velocities  are  Vj ,  Vj2 ,  Vjj ,  and  the  z  velocities  are 
w. ,  Wj2 ,  Wi3 .  For  layers  2  and  3,  the  x,  y,  and  z  velocities  at  node  i  are 


Ui2=u,+H0- 

(20) 

Ui3=Ui-H0l 

(21) 

Vi2=Vi 

(22) 

Vi3  =  Vi 

(23) 

Wi2  =Wi-H0; 

(24) 

Wi3=Wi-HH0‘ 

(25) 

The  offset  between  layers  is  simply 

H  =  V„/3A  (26) 

where  Vo  is  the  initial  volume  and  A  is  the  current  area  of  the  triangular  element.  The  x,  y,  and  z 
velocities  for  nodes  j  and  k  can  be  obtained  in  a  similar  manner. 

The  in-plane  strain  rates  (e„,  e^,  y^z)  obtained  in  exactly  the  same  manner  as  for  2D 
triangular  elements.  See  Equations  1 1, 12,  and  14  in  subsection  3.4.4.  The  rotational  velocities 
are  assumed  to  vary  hnearly  between  the  nodes  in  the  same  manner  as  the  translational  velocities. 
The  final  strain  rates,  for  the  local  coordinate  system  in  Figure  21,  are  as  follows: 
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K-- 

Xk-Xj)(wi-y0O-Xk(wj 

(28) 

=  f(e.. 

^y) 

(29) 

1  I 

/  \1  1/Ai  Ai  A 

Yxy 

■2aI 

.Zk(vj-Vi)J-l-j(0‘ -J- 0^,-1- 0 

z) 

(30) 

Yxz 

(x,-x,)v,-x^v,  +  XjV,]- 

j(e’x+ei+0!:) 

(31) 

-y0x)~^k(wi-y0O+(xk 

-x.)(uj-t-y0‘)-Xk(uj+y0i)  +  Xj 

(u,-hy0^)](32) 

The  through  thickness  strain  rate,  £y ,  cannot  be  determined  explicitly,  but  must  be  iterated  to 

provide  <Ty  =  0  through  the  thickness.  Also  note  that  for  layer  1  (center)  y  =  0,  for  layer  2  (top)  y 
=  H,  and  for  layer  3  (bottom)  y  =  -H. 

The  determination  of  the  in-plane  shear  and  deviator  stresses  is  similar  to  that  used  for  the  2D 
bending  shells.  Because  the  stress  normal  to  the  plane  is  Cy  =  0,  an  iteration  procedure  must  be 
used.  The  iteration  begins  by  assuming  incompressible  flow,  which  gives 

ey=-(e.+ez)  (33) 


Now  that  all  the  strain  rates  are  defined  (ex>^y>^z»Yxy>Yxz»yyz)  shear  and  deviator  stresses 

can  be  determined  in  the  standard  manner,  as  provided  in  subsection  3.8.4.  The  iteration  on  the 

volumetric  strain,  to  obtain  ay  =  0,  is  almost  identical  to  that  presented  in  subsection  3.4.6  and 
shown  in  Figure  12.  For  2D  bending  shells  the  known  normal  strain  rates  are  e'  and  ,  and  the 

through  thickness  strain  rate  (to  be  determined)  is  e' .  For  3D  bending  shells  (as  well  as  plane 
stress  elements),  the  known  normal  strain  rates  are  and  e^,  and  the  through-thickness  strain 
rate  (to  be  determined)  is  .  Other  than  the  change  in  nomenclature,  the  algorithms  for  the 

three  cases  (2D  bending  shells,  plane  stress  elements,  3D  bending  shells)  are  identical. 

The  stresses  acting  on  the  elements  are  shown  in  Figure  24,  and  the  resulting  forces  are  shown  in 
Figure  25.  The  x,  y,  and  z  forces  on  nodes  i,  j,  and  k  (for  each  layer)  are  as  follows: 
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F'=f[T„(x,-xJ-a,z,] 

(34) 

F;=f[vK-Xj)-v,] 

(35) 

Fz  =  Y[<yz(xk  “Xj^-T^j^Zk] 

(36) 

i  Hr  1 

(37) 

Fy=f  [Vk-Vk] 

(38) 

i  Hr  1 

Fi=YKZk-<yzXkJ 

(39) 

Fx  =f  (Vj) 

(40) 

il 

(41) 

(42) 

where  the  stresses  are  for  the  specific  layer  being  considered. 

Because  there  are  shear  forces  in  each  layer,  it  is  necessary  to  have  concentrated  moments  on  the 
nodes  of  each  layer  to  provide  equilibrium.  Summing  moments  about  the  z  axes  at  node  i  gives 

3AM,  +  Fy"x,+Fjxj=0  (43) 

The  resulting  concentrated  moments  are 

AM^  =  -(Fy  Xk  +  Fy^Xj)  /  3  (44) 
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Simiming  about  the  x  axis  at  node  i  gives 

AM,=FX/3  (45) 

Again,  the  concentrated  moments  for  each  layer  must  use  the  appropriate  force  for  the  specific 
layer  being  considered. 

The  next  step  is  to  sum  the  forces  and  moments  firom  each  layer  onto  the  three  nodes. 

For  node  i,  the  three  forces  are 


-  Px{l)'^Px(2)‘^^x(3) 

(46) 

Py  =Fy(0  +  Py(2)  +  Fy(3) 

(47) 

Fz=F-,)  +  F',)  +  F‘3) 

(48) 

In  Equation  34,  F^  is  for  a  generalized  single  layer,  but  in  Equation  46  it  represents  the  forces  for 
all  three  layers.  F^(,)  in  Equation  46  is  F^  from  Equation  34  for  layer  1.  The  forces  on  nodes  j 

and  k  are  determined  in  a  similar  manner. 

The  net  moments  on  node  i  are 

M'  =  AM.,,,  +  AM.„,  +  AM.O,  +  F;„,H  -  f;,„H  (49) 

M'  =  am.,„+am„„+am.,„-f;,„h+f;,„h  (50) 

The  net  moments  on  nodes  j  and  k  are  determined  in  a  similar  maimer. 

The  final  step  is  to  convert  the  forces  and  moments  back  to  the  system  coordinates,  from  the 
local  coordinate  system. 

The  heat  conduction  algorithm  is  identical  to  the  2D  heat  conduction  algorithm  for  triangular 
elements.  For  the  heat  conduction  option  all  three  layers  will  have  equal  temperatures  because 
all  layers  are  attached  to  the  same  nodes. 
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3.8.8  Membrane  Shell  Elements 

The  3D  membrane  shell  elements  are  a  simplified  form  of  the  3D  bending  shell  elements 
presented  in  subsection  3.8.7  and  shown  in  Figure  21.  The  membrane  shell  element  has  only  in¬ 
plane  stresses  and  is  composed  of  only  a  single  layer.  It  cannot  produce  bending  or  transverse 
shear  stresses.  Furthermore  it  is  assumed  to  be  incompressible,  and  this  eliminates  the  need  to 
iterate  to  obtain  the  through-thickness  stress.  This  algorithm  is  presented  in  Reference  21  for 
plane  stress  geometry. 

Referring  back  to  Equation  29  of  subsection  3.8.7,  the  through-thickness  strain  rate  for 
incompressible  deformation  is 


where  and  are  the  in-plane  strain  rates  in  the  local  coordinate  system. 

The  deviator  stresses  (s,^,  Sy ,  s^j  can  then  be  determined  from  e^,  Ey,  6^ ,  which  gives  the  net 
stresses 


a,=s,-P  (2) 

Oy  =  Sy  -  P  (3) 

a,  =  s,-P  (4) 

In  Equation  3,  setting  ay  =  0  gives  P  =  Sy,  and  this  provides  an  explicit  determination  of  the 
stresses  in  Equations  2—4. 

3.8.9  Nonreflective  Boundary  Elements 

The  3D  nonreflective  boundary  elements  are  incorporated  as  infinitely  thin,  triangular,  massless 
elements  as  shown  in  Figure  26.  They  can  be  used  to  absorb  elastic  waves  to  represent  infinite 
media,  and  they  can  also  be  used  to  reduce  the  size  of  the  model  for  other  applications  such  as  for 
penetration  problems  (References  8-10). 
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Figure  26.  Description  of  the  3D  Nonreflective  Boundary  Eiement 

The  average  velocity  normal  to  the  triangular  face  is 

(1) 

where  u  =  +Uj  +u,j)/3  is  the  average  velocity  in  the  x  direction,  and  v  and  Iv  are  the 
average  velocities  in  the  y  and  z  directions.  The  direction  cosines  normal  to  the  surface  are 


f,=A,/A 


(2) 


iy=Ay/A 


(3) 
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£,  =  AJA 


(4) 


where  Ax,  Ay,  and  Az  are  the  projected  areas  on  the  yz,  xz,  and  xy  planes  when  viewed  along  the 
X,  y,  and  z  axes,  and  A  is  the  total  area. 


A.  41 

(yi-yj)(^k-Zi)-(yk-yj)(2i-^j)] 

(5) 

> 

II 

ts>| 

(6) 

1 

^-2 

(x,  -  Xj)(yk  -  yi)-(xk  -  Xj)(y,  -  yj)] 

(7) 

a=7a 

(8) 

In  addition  to  the  normal  velocity  it  is  necessary  to  determine  two  velocities  in  the  plane  of  the 
element.  The  first  of  the  in-plane  velocities  is  in  the  direction  along  the  line  connecting  nodes  j 
and  k. 


V,=U^„+V^y,+Wf„ 

where  Vi  is  the  in-plane  velocity,  h,  V,  and  W  are  the  average  nodal  velocities  as  defined 
previously,  and  the  directions  cosines  are 


^xi=(xk-Xj)/d 

(10) 

(11) 

(12) 

where  d  is  the  distance  between  nodes  j  and  k. 

The  other  in-plane  velocity  is  perpendicular  to  both  Vn  and  Vi.  It  is  expressed  as 
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(13) 


^2  =U^x2  +  V^y2+W^z2 

where  the  direction  cosines  are 
/  =P  f  -P  P 

^  x2  ^  zl  ^  z^  yl 

P  zzP  P  -P  P 
^  y2  ^  z^  xl  ^  x^  zl 


The  total  normal  resisting  force  on  the  element  is 


(14) 

(15) 

(16) 

(17) 


and  the  total  in-plane  forces  are 

^1  ~  ~Po^shear-'^% 

(18) 

F2  =-poCshearAV2 

(19) 

where  po  is  the  initial  density  of  the  material,  A  is  the  total  area  and  Vn,  Vi  and  N 2  are  the 
velocities.  The  longitudinal  and  shear  soimd  velocities  are 


c.  =  .J(K,+4G/3)/p. 

(20) 

^shear 

(21) 

where  Ki  is  the  bulk  modulus  and  G  is  the  shear  modulus. 

Now  the  total  forces  must  be  aligned  with  the  three  principal  axes  and  distributed  equally  to  the 
three  nodes. 

F;  =  Fi  =  F;  =  (f^4  +  F/.,  +F,£.,)/3  (22) 
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F;=F.j  =  Fy  =(Vy  +  F.^y.+F2^y2)/3 

F‘  =Fj  =  =(Fn4  +  F/z.  +F2^z2)/3 


(23) 

(24) 
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SECTION  4 
SPH  ALGORITHMS 


During  the  past  few  years  SPH  (Smooth  Particle  Hydrodynamics)  methods  have  been  developed 
and  applied  to  problems  involving  high  velocity  impact.  The  appeal  of  SPH  for  high  velocity 
impact  is  that  it  is  a  Lagrangian  technique  and  that  it  has  variable  connectivity  to  allow  for  severe 
distortions.  The  Lagrangian  feature  is  desired  because  it  allows  the  grid  to  be  embedded  in  the 
material  and  this  reduces  some  of  the  material  interface  and  material  history  problems  associated 
with  Eulerian  codes.  Furthermore,  the  ability  to  handle  severe  distortions  allows  the  SPH 
technique  to  be  applied  to  problems  that  historically  have  been  reserved  for  Eulerian  approaches. 

Although  SPH  approaches  can  be  applied  to  severe  distortions,  they  are  generally  not  as  good  as 
standard  finite  elements  for  structural  response  applications.  These  same  characteristics  are  held 
by  Eulerian  techniques.  The  important  advantage  of  SPH,  however,  is  that  its  Lagrangian 
formulation  allows  it  to  be  straightforwardly  linked  to  standard  finite  element  Lagrangian 
formulations.  This  means  that  it  is  possible  for  both  severe  distortions  and  structural  response 
computations  to  be  performed  with  a  single  Lagrangian  code,  and  that  both  severe  distortions  and 
structural  responses  can  occur  in  the  same  problem.  A  long  term  objective  is  to  allow  the  user  to 
define  almost  any  impact  problem  with  a  standard  finite  element  grid,  and  then  to  allow  the 
standard  elements  to  be  converted  to  SPH  nodes  as  the  standard  elements  become  distorted. 
Although  this  approach  has  been  demonstrated,  more  work  is  required  to  increase  the  accuracy 
and  robustness  for  a  wider  range  of  problems.  The  algorithms  that  follow  are  reported  in 
References  28-31. 

4.1  2D  AXISYMMETRIC  GEOMETRY 

A  schematic  overview  of  the  structure  of  the  EPIC  code  is  shown  in  Figure  27.  It  is  very  similar 
for  both  standard  elements  and  SPH  nodes,  with  the  primary  differences  being  the  computations 
of  the  strains,  strain  rates  and  nodal  forces.  The  determination  of  the  nodal  displacements  and 
velocities,  as  well  as  the  stresses,  is  identical  for  both  approaches  and  is  not  included  here.  A 
description  of  the  basic  2D  axisymmetric  finite  element  algorithm  is  provided  in  subsection  3.4. 

Figure  28  represents  some  features  of  the  SPH  technique.  Node  i  is  designated  as  the  center  node 
and  the  neighbor  nodes  are  designated  as  nodes  j.  The  distance  between  nodes  is  ry,  the 


C  t19616.doc 


117 


Current  Velocities  and 
Displacements  of  all  Nodes 


Standard  Elements 

Determine  Strains,  Strain 
Rates,  etc.,  in  all  Elements 


SPH  Nodes 

Determine  Strains,  Strain 
Rates,  etc.,  In  all  SPH  Nodes 


Material  Models 


Determine  Stresses  from 
Strains,  Strain  Rates,  etc. 


diameters  of  the  nodes  are  di  and  dj,  and  the  masses  of  the  nodes  are  Mi  and  Mj.  The  masses 
remain  constant  throughout  the  computation,  and  are  obtained  from  M  =  poVo  where  po  and  Vo 
represent  the  initial  density  of  the  material  and  the  initial  volume  represented  by  the  node. 

The  SPH  approach  allows  for  variable  nodal  connectivity,  and  this  means  that  it  is  necessary  for 
each  center  node  i  to  search  and  find  the  closest  neighbor  nodes  j.  This  searching  time  can  be 
significant,  especially  for  larger  problems,  and  several  approaches  have  been  tried.  The  EPIC 
code  currently  uses  a  bucket  searching  algorithm  (Reference  13). 

4.1.1  Smoothing  Functions 

The  smoothing  function  is  an  important  part  of  the  SPH  algorithm.  Two  smoothing  functions 
and  their  derivatives  are  shown  in  Figure  29.  These  smoothing  functions  can  be  used  for  both 
plane  strain  and  axisymmetric  geometries.  The  B-Spline  has  been  a  commonly  used  smoothing 
function  in  the  past,  but  the  recently  introduced  Quadratic  smoothing  function  appears  to  have 
some  advantages  (Reference  30). 

The  B-Spline  smoothing  ftmction  is  expressed  as 

where  uy  =  ry/hy,  and  the  smoothing  distance  is 

h,  =  a(d,  +  d,)/2  (2) 

The  dimensionless  smoothing  distance,  a,  is  a  user  supplied  input.  It  is  usually  taken  as  a  = 

1 .0,  but  other  values  can  be  used  (0.8  <  a  <  1 .2).  The  diameters,  di  and  dj,  can  be  obtained  (for 
axisymmetric  geometry)  from 


r.v 
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Dimensionless  Distance  Between  Nodes,  Vjj  =  ry/  hy 
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Figure  29.  2D  B-Spline  and  Quadratic  Smoothing  Functions  and  Their  Derivatives 
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where  8v  is  the  volumetric  strain  (defined  later),  do  is  the  initial  node  diameter,  and  Xo  and  x  are 
the  initial  and  current  x  (radial)  coordinates. 


A  requirement  of  the  smoothing  fimction  is  that  it  exhibits  the  characteristics  of  a  Dirac  delta 
fimction  as  hij  approaches  zero.  It  will  be  shown  later  that  it  is  the  derivative,  ,  of  the 
smoothing  function  that  actually  provides  the  weighting  function  for  the  strain  rates  and  forces. 
The  negative  value  of  the  B-Spline  derivative,  -  ,  is  shown  in  Figure  29  and  is  expressed  as 


TthJ 


30  45 

- 1)  H - 

7  «  14 


t)' 

ij 


0  <  Uij  <  1 


(4a) 


1  <  Uij  <  2 


(4b) 


The  interesting  feature  of  the  B-Spline  derivative  is  that  -  Wj-  exhibits  a  maximum  at  oy  =  2/3 
and  decreases  for  both  smaller  and  larger  values  of  oy.  For  uy  >  2/3  it  is  logical  that 
decreases  because  it  has  less  influence  as  the  distance  is  increased.  To  have  -  W/  decrease  for 
Uy  <  2/3  is  not  intuitively  satisfying  because  the  closer  neighbor  node  j  comes  to  center  node  i, 
the  less  influence  it  has.  This  can  also  lead  to  instabilities  in  compression  (Reference  32). 


The  Quadratic  smoothing  function  is  also  shown  in  Figure  29.  It  is  expressed  as 


8^‘j  2^‘-''*'2 


0  <  Uy  <  2 


(5) 


The  derivative  is 


3 

2 


0  <  Uy  <  2 


(6) 


For  the  Quadratic  smoothing  function  derivative,  the  weighting  function,  ,  always  increases 

as  the  nodes  move  closer  together,  and  always  decreases  as  they  move  apart.  This  intuitively 
appears  to  be  more  realistic  than  the  B-Spline  derivative.  It  is  also  more  simple  than  the  B- 
Spline  and  does  not  have  the  compressive  softening  that  can  lead  to  compressive  instabilities. 


Based  on  these  arguments  alone,  it  would  appear  that  the  newer  Quadratic  smoothing  function  is 
an  improvement  over  the  traditional  B-Spline. 

4.1.2  Strain  Rates 

For  2D  axisymmetric  geometry,  the  three  normal  strain  rates  ,  6^ ,  £9 ) ,  the  shear  strain  rate, 

,  the  rotational  rate,  coxz,  and  the  volumetric  strain  rate,  ,  for  center  node  i,  are  as  follows. 


£.=-Z|3.W^  V,(Uj-a,K,/2ia, 

j 

(1) 

£.=-XP,Wi' V,(Vj-ViK./2ltx, 

j 

(2) 

(3) 

=  -I  -  il, )« .  +  P.(vj  -  ]  /  2icxj 

j 

(4) 

<i>„  =Sw^j[p.(uj -u,)f , -v,>,]/4)tx, 

j 

(5) 

e  =£,+£,+£9 

(6) 

where  =  9Wy  /  9r  is  the  derivative  of  the  smoothing  function,  Vj  is  the  current  volxime  of  the 
nodej,  Uj  and  Uj  are  the  x  velocities  ofnodesi  and  j,  Vj  and  Vj  are  the  z  velocities,  and 

are  the  direction  cosines  from  node  i  to  nodej,  and  xj  is  the  x  coordinate  of  node  j.  The  updated 
volumetric  strain  is  obtained  by  integrating  the  volumetric  strain  rate. 

£v‘''^‘ =£l+elAt(l+£0 

where  At  is  the  integration  time  increment,  and  the  factor  (l  +  £t )  converts  the  strain  rate  from 
the  current  configuration  back  to  the  initial  configuration. 

The  three  p  factors  are  used  to  normalize  the  smoothing  functions  such  that  they  will  provide  the 
exact  strain  rates  in  the  three  principal  directions  for  states  of  constant  strain  rates.  They  are 
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obtained  from  Equations  1-3  by  setting  Uj  -  lij  =  ^  constant  strain  rate  in  the  x 

direction,  Vj  -  Vj  s  e^^Tij^z  ®  constant  strain  rate  in  the  z  direction,  and  Uj  =  e^Xj  for  a  constant 

hoop  strain  rate.  The  resulting  normalizing  factors  are 


Sw;Vjr/,/2jlXj 

_ -1 


(8) 

(9) 

(10) 


The  effect  of  using  the  Normalized  Smoothing  Function  (NSF)  algorithm  can  be  significant. 
Figure  30  shows  a  cross-section  of  an  axisymmetric  ring  of  material  represented  by  25  SPH 
nodes  (5x5).  The  two  cross-sections  on  the  left  side  are  for  the  Standard  Smoothing  Fimctions, 
which  means  that  they  do  not  use  the  NSF  algorithm.  The  two  on  the  right  side  do  use  the  NSF 
algorithm.  The  top  two  use  the  B-Spline  smoothing  function  and  the  bottom  two  use  the 
Quadratic  smoothing  function.  The  numbers  in  the  centers  of  the  circular  SPH  nodes  represent 
the  equivalent  strain  rate  in  the  node  when  the  axisymmetric  cross-section  is  subjected  to  a 
stretching  radial  velocity  of 


Uj  =  3Xj  /  2 


(11) 


where  xj  is  the  radial  coordinate.  The  resulting  uniform  strain  rates  are 


=  du/8x  =  3/2 

(12) 

=  u/x  =  3/2 

(13) 

o 

II 

II 

(14) 
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Figure  30.  Equivalent  Strain  Rates  in  Radially  Stretching  Cross-Sections  With  a 
Uniform  Arrangement  of  SPH  Nodes 


The  equivalent  strain  rate  is  defined  as 
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Substituting  Equations  12-14  into  Equation  15  gives 

is  1.0  (16) 

The  upper  left  cross-section  in  Figure  30  (B-Spline  without  NSF)  shows  good  acciiracy  in  the 
center  region  (e  =  1.0 1) ,  but  the  side  boundaries  and  comers  are  significantly  low 

(0.43  <  8  <  0.85) .  The  reason  for  the  low  strain  rates  at  the  boimdaries  is  that  the  basic  SPH 

algorithm  is  for  an  interior  node,  which  assumes  there  is  a  full  distribution  of  neighbor  nodes. 
When  the  distribution  of  neighbor  nodes  is  not  full,  inaccuracies  are  introduced.  The  lower  left 
in  Figure  30  (Quadratic  without  NSF)  has  a  similar  pattern,  but  the  interior  nodes  have  a  lower 
strain  rate  (e  <  0.8?) .  This  does  not  necessarily  mean  that  the  B-Spline  is  more  accurate  than  the 

Quadratic  smoothing  function,  because  the  results  in  Figure  30  are  specifically  for  a  uniform  grid 
with  a  dimensionless  smoothing  distance  of  a  =  1.0.  For  some  other  conditions  (such  as  a  =  0.8) 
the  Quadratic  smoothing  function  is  more  accurate,  as  will  be  shown  later.  At  larger  smoothing 
distances  the  interior  nodes  approach  an  equivalent  strain  rate  of  e  =  1.0 ,  and  at  lower  smoothing 
distances  (a  <  0.8)  the  strain  rates  decrease  to  8  =  0  at  a  =  0.5.  This  effect  is  illustrated  in 
Figme  31,  where  the  equivalent  strain  rate  of  an  interior  node  is  shown  as  a  function  of  the 
dimensionless  smoothing  distance,  a.  Note  that  the  B-Spline  is  more  accurate  at  a  =  1.0  and  the 
Quadratic  is  more  accurate  at  a  =  0.8.  The  sharp  slope  changes  for  the  Quadratic  smoothing 

function  at  a  =  0.707,  1.000, 1.118,  etc.,  occur  when  the  expanding  smoothing  distance 
encoimters  additional  nodes  at  Ojj  =  2.0.  The  Quadratic  smoothing  function  derivative,  W^' ,  has  a 

finite  slope  at  oy  =  2.0  (and  therefore  has  a  noticeable  effect)  whereas  the  slope  of  for  the  B- 

Spline  is  zero  at  uy  =  2.0.  When  the  NSF  algorithm  is  applied,  as  shown  on  the  right  side  of 
Figure  30  and  in  Figure  31,  then  the  results  are  much  improved.  The  NSF  algorithm  is  also  very 
effective  when  the  nodes  are  in  a  non-iuiiform  arrangement,  and  it  is  shown  to  provide  excellent 
results  for  several  cylinder  impact  examples  (Reference  30). 

4.1.3  Forces 

After  the  strain  rates,  rotational  rate,  and  volumetric  strain  in  Equations  1-7  of  subsection  4.1.2 
are  determined,  it  is  possible  to  determine  the  shear  and  deviator  stresses,  the  pressure,  and  the 
nodal  artificial  viscosity  in  the  standard  manner,  as  provided  in  subsection  3.4.4.  These  stresses 
must  then  be  converted  to  forces  as  shown  in  Figure  27. 
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Figure  31 .  Equivalent  Strain  Rate  Versus  Dimensionless  Smoothing  Distance  for 
an  interior  SPH  Node  in  a  Radial  Stretching  Cross-Section 


The  nodal  force  in  the  x  direction  on  node  j,  due  to  the  stress  of  node  i,  is 


pij  —  py  (P^ane)  ^  pU  (hoop) 


(1) 


The  force  due  to  the  in-plane  stresses  is 

p«(p.ane)  =  Wi;Vi  Vj[p,  (<  -  Qy)f ,  +  PXf,]/ (2) 

where  =  sj^  -  (Pj  +  Qj )  is  the  net  normal  stress  in  the  x  direction,  composed  of  the  deviator 
stress,  pressure,  and  nodal  artificial  viscosity,  and  is  the  shear  stress. 

There  is  also  an  artificial  viscosity,  Qy,  which  is  dependent  on  the  relative  velocities  of  nodes  i 
and  j.  It  is  intended  to  stabilize  the  grid  and  keep  adjacent  nodes  fi-om  becoming  too  close  to  one 
another.  This  will  be  designated  as  a  bond  viscosity,  because  it  acts  on  the  bond  between  nodes  i 
and  j.  It  can  also  introduce  a  significant  amount  of  artificial  strength. 
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The  force  due  to  the  hoop  stress  is 


pu  (hoop)  ^  /  47tr/  (3) 

and  the  force  in  the  z  direction  is 

pj  ='W|;V,V,[p.(<-Q,K+P.t|„<.]'2)tXj  (4) 

Equations  2-4  provide  the  forces  only  on  the  neighbor  nodes  j.  The  forces  on  center  node  i,  due 
to  the  stresses  in  node  i,  are  equal  and  opposite  to  the  in-plane  forces  in  the  x  and  z  directions. 


E 

i 

pij  (plane) 

(5) 

I 

(6) 

j 


4.1.4  Artificial  Viscosity 

There  are  two  forms  of  artificial  viscosity;  the  nodal  viscosity,  Qi,  and  the  bond  viscosity,  Qy. 

The  important  distinction  between  Qi  (nodal)  and  Qy  (bond)  is  that  Qi  can  only  be  activated 
during  volumetric  strain  rate,  but  Qy  can  be  activated  even  if  there  is  no  volumetric  strain  rate 
(such  as  pure  shear  or  incompressible  flow).  The  net  effect  is  that  Qy  can  introduce  additional 
(artificial)  strength  into  the  computed  results. 

The  nodal  viscosity  is  identical  to  that  used  in  standard  finite  element  and  finite  difference 
methods  (References  5  and  33): 

Q,=CtPiC|h||E.|+CgP|hf£;  (1) 

for  <  0.  Cl  and  Cq  are  the  linear  and  quadratic  coefficients,  pi  is  the  density  of  node  i,  Ci  is  the 
sound  velocity  of  node  i,  and  hi  is  the  minimum  smoothing  distance  between  node  i  and  neighbor 
nodes). 
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The  bond  viscosity  of  Monaghan  and  Gingold  (Reference  34)  can  be  expressed  in  the  following 
form: 


Qu  =  CLpiCi|Hy|  +  CQPiH? 


for  pij  <  0,  and 


=  _2 


h.j(V.j  +  Vij) 


lij+eohfj 


(2) 


(3) 


where  Uj  =  U;  -  Uj  is  the  x  velocity  difference  between  nodes  i  and  j,  and  xy  =  xj  -xj  is  the  x 
coordinate  difference.  The  analogous  terms  for  the  z  components  are  Vjj  and  zy.  The  distance 

between  nodes  is  rjj  =  +  Zy ,  and  So  is  a  small  number  (So  «  0.01)  that  acts  to  limit  py  as  ry 

becomes  small. 

Equation  2  can  be  rewritten  in  the  following  form  (for  8o  =  0): 

Qij  =  CLpiCihjjEiJ  +  CQPihfjefj  (4) 

for  £;•  <  0.  This  looks  very  similar  to  the  nodal  artificial  viscosity  of  Equation  1,  except  the 
volumetric  strain  rate,  e„ ,  of  Equation  1  is  replaced  with  a  linear  strain  rate,  Ey ,  in  Equation  4. 

This  linear  strain  rate  is  simply  the  strain  rate  along  the  bond  between  nodes: 

where  is  the  velocity  of  node  i  along  the  bond  from  node  i  to  node  j,  and  Vf  is  the  velocity 

of  node  j  in  the  same  direction.  A  similar  form  of  this  bond  viscosity  was  developed  previously 
for  the  NABOR  particle  method  algorithm  (References  35  and  36). 

It  should  be  noted  that  most  of  the  SPH  literature  (Reference  34)  describes  the  bond  viscosity 
term  as 


(6) 
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for  (iy  <  0.  Here  Cl,  Cq,  and  ^y  are  as  defined  previously,  cy  is  the  average  sound  velocity  at 
nodes  i  and  j,  and  py  is  the  average  density  at  nodes  i  and  j.  Also,  Oy  is  the  viscosity  for  both 
nodes  i  and  j.  In  the  algorithms  of  Equations  1-5,  Qy  is  based  on  the  density  and  sound  velocity 
of  node  i  only.  The  relationship  between  the  two  forms  is  approximated  by 

Q«"ni|p?/2  (7) 

For  identical  sound  velocities  and  densities  at  nodes  i  and  j,  the  relationship  of  Equation  7  is 
exact.  The  factor  2  in  the  denominator  occurs  because  fly  is  for  both  nodes  whereas  Qy  is  only 
for  node  i  (and  there  is  another  Qji  at  node  j). 

4.1.5  Material  Interfaces 

Figure  32  shows  a  material  interface  represented  by  SPH  nodes.  If  the  two  materials  are  not 
bonded  together,  the  standard  SPH  algorithm  introduces  unacceptably  large  errors  because  nodes 
from  material  A  influence  the  strain  rates  in  nodes  from  material  B,  and  vice  versa.  Also,  shear 
and  tensile  stresses  are  developed  between  the  two  materials  such  that  sliding  and  separation  are 
significantly  inhibited.  If  SPH  approaches  are  to  be  applied  to  problems  involving  sliding 
interfaces,  then  specialized  interface  algorithms  must  be  developed. 


>  Materia!  A 


Figure  32.  Material  interface  With  SPH  Nodes 


4.2  2D  PLANE  STRAIN  GEOMETRY 


The  2D  plane  strain  geometry  is  a  simplified  form  of  the  2D  axisymmetric  geometry.  The 
primary  difference  is  that  the  plane  strain  geometry  is  for  a  unit  thickness  and  has  no  hoop  strains 
or  strain  rates. 

Using  the  same  notation  as  for  the  axisymmetric  geometry  in  subsection  4.1.2,  the  two  normal 
strain  rates  (e^,  e^),  the  shear  strain  rate,  ,  and  the  rotation  rate,  coxz,  for  center  node  i,  are  as 

follows; 


K  =-SP.W^i(4j-4iK 

j 

(1) 

e,  =  -SP.W^j(vj-v,)4 

j 

(2) 

Yxz  = -E  WifVj  [pz  (u j  -  Ui  K  +  Px  (v j  -  Vi  ] 

j 

(3) 

=EW/V^[p.(v.,  -f,)^]/2 

(4) 

J 


One  of  the  differences  is  that  the  volume  of  node  j  is  now  identical  to  the  area  (Vj  =  Aj)  because 
the  plane  strain  geometry  is  for  a  unit  thickness.  Also,  the  axisymmetric  hoop  strain  rate  (ej)  is 

replaced  by  the  through-thickness  strain  rate,  which  is 

(5) 


The  normalizing  factors  are  also  modified. 
a  -1 

R  -1 


(6) 

(7) 
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After  the  stresses  have  been  determined,  the  forces  on  node  j,  due  to  stresses  on  node  i,  are  as 
follows: 


P#  =  W^V,Vj[|3.(<  -Q,)« . 

P»  =w;ViVj[p.(ai-QsK  +  PX4] 

Here  the  notation  is  as  provided  in  subsection  4.1.3. 


4.3  3D  GEOMETRY 


There  is  very  little  additional  complexity  in  going  from  2D  to  3D  geometry.  Again,  there  is 
similarity  between  the  3D  SPH  and  the  3D  standard  element  algorithms  provided  in  subsection 
3.8.  The  B-Spline  smoothing  function  in  3D  is 


0  <  Uij  <  1 


(la) 


1  <  Uij  <  2 


(lb) 


and  the  derivative  is 


0  <  Uij  <  1 


(2a) 


W;'  = 


1 


1  <  Uij  <  2 


(2b) 


The  3D  B-Spline  smoothing  function  and  its  derivative  are  simply  7/lOhij  times  the 
corresponding  2D  functions  given  in  Equations  1  and  4  of  subsection  4.1.1,  and  shown  in 
Figure  29. 

The  3D  form  of  the  Quadratic  smoothing  function  is 
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0  <  Oij  <2 


(3) 


5(3  2 


3  '  3) 


and  the  derivative  is 


W'  =  — — 

_u4 


1  ri5fi  , 

■  —  V  -  -1 

;th’Ll6l2 


O^-Uij  <2 


(4) 


Here  again,  the  3D  Quadratic  smoothing  function  and  its  derivative  are  simply  S/Shy  times  the 
corresponding  2D  functions  given  in  Equations  5  and  6  of  subsection  4.1.1,  and  shown  in 
Figure  29. 

The  three  normal  strain  rates  (e^  ,ey  ,e  J ,  the  three  shear  strain  rates  (Y,^y  ,y^ , Yy^ ) ,  and  the  three 
rotational  rates  (©,jy  ,co,(^  ,©y^ )  for  the  center  node  i,  are  as  follows: 


e,=-SP,W,;  V/Uj-u,K, 

J 

j 

8,  =-IP,Wj'  Vj(Wj-WiK, 

j 

■y„  =-ZW|' Vj[p,(u^-u,K,+p.(Vi-v,)<j 

Y„  =-IW/  Vj[P.(Wj  -w,K,  +p.(u,-u,K.] 

j 

=  -I  W,'  Vj[p. (Vj  -  V,  K ,  +  P, (Wj  -  w, )^ , ] 

0),,  =ZW|'  Vj[P,(u,-u,K,-P,(Vj-V|KJ/2 

<»„  =  E  W;  Vj[P,(Wj  -  w,)<.  -P.(il,  -  u,K.]/2 

j 

(Oy,  =  iw^  Vj[p,(v.  -  Vi)^,  -Py(w.  -  WiKy]/2 

j 


(5) 

(6) 

(7) 

(8) 
(9) 

(10) 

(11) 

(12) 

(13) 
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There  is  a  change  in  nomenclature  for  the  3D  algorithm  when  compared  to  the  2D  algorithm.  For 
the  2D  strain  rates  and  rotational  rates  in  Equations  1-5  of  subsection  4.1.2,  U;  and  Vj  are  the 

velocities  corresponding  to  the  x  (radial)  and  z  (axial)  directions,  respectively.  For  the  3D 
algorithm  of  Equations  5-13,  the  x,  y,  and  z  velocities  are  represented  by  U; ,  Vj ,  and  W; , 

respectively.  This  is  consistent  with  the  3D  finite  element  algorithms  of  subsection  3.8. 


The  three  p  factors  are  obtained  in  the  same  manner  as  the  2D  geometry  and  are  expressed  as 


Px  = 

Py  = 


-1 


-1 

-1 

SK  V.r/, 


(14) 


(15) 

(16) 


There  is  an  additional  complexity  in  applying  the  NSF  algorithm  in  3D  geometry.  For  2D 
axisymmetric  geometry  there  is  only  one  possible  orientation  of  axes  because  the  z  axis 
represents  the  axis  of  rotation.  For  3D  geometry  the  p  factors  are  determined  along  the  three 
principal  axes.  If  the  axes  are  changed,  then  the  p  factors  may  change  and  the  computed  results 
may  also  change.  Additional  work  is  required  to  examine  the  magnitude  of  this  potential 
problem. 

Finally,  the  nodal  forces  are  as  follows: 


Pi* -Qij)4 +Pz<^z] 
Pi*  =  W*;  XVj[P,(a;  -Q,)^,  +PX/x  +Pz<4] 
Pi*  =  K  V<VJP,(al  -Q,)^,  +PX4  +P,x;/J 


(17) 

(18) 
(19) 
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SECTION  5 

UNKING  OF  FINITE  ELEMENTS  AND  SPH  NODES 


It  has  been  noted  previously  that  a  desirable  characteristic  of  the  SPH  approach  is  that  it  can  be 
linked  to  standard  finite  element  grids.  This  allows  highly  distorted  materials  to  be  considered 
together  with  structural  response  materials.  Also,  as  shown  in  Figure  27,  a  single  material  model 
subroutine  can  be  used  for  both  the  standard  elements  and  the  SPH  nodes.  This  means  that  a 
material  model  needs  to  be  incorporated  and  validated  only  one  time,  and  that  the  user  can  be 
assured  that  the  identical  material  model  is  used  for  both  the  standard  elements  and  the  SPH 
nodes.  This  work  has  been  reported  in  References  28, 29,  and  3 1 . 

Figure  33  shows  how  SPH  nodes  can  be  attached  to  a  standard  finite  element  grid.  The  broken 
circle  around  node  i  represents  the  effective  region  (0  <  oy  <  2)  for  SPH  nodes  on  node  i.  The 
size  of  the  SPH  interface  nodes  is  equal  to  the  size  of  the  interface  elements.  The  mass  of  the 
SPH  interface  nodes  comes  from  the  SPH  nodes  only,  and  this  ensures  that  the  mass  on  all  nodes 
is  essentially  equal. 

The  current  approach  for  this  type  of  linkage  is  to  determine  the  strain  rates  for  SPH  node  i  only 
fi-om  nodes  ni, ...,  ns.  The  standard  nodes  in  the  effective  region  (ns, ...,  ng)  are  not  included 
because  of  the  complexity  they  would  introduce  into  the  numerical  algorithm.  The  strain  rates 
can  be  accurately  computed  on  the  interface  by  using  the  NSF  algorithm,  but  the  resulting  forces 
appear  to  be  in  error. 

The  forces  on  SPH  nodes  i  come  fi-om  SPH  nodes  ni, . . .,  ns  and  fi-om  interface  elements  B  and 
C.  There  is  no  direct  force  contribution  fi-om  standard  nodes  ns, . . .,  n9,  other  than  through  the 
interface  elements.  Improvements  could  include  consideration  of  standard  nodes  (ns, . . .,  ng)  or 
generation  of  ghost  nodes  on  the  standard  node  side  of  the  interface. 

Figure  34  shows  how  SPH  nodes  interact  with  a  standard  grid  on  a  sliding  interface.  The  SPH 
nodes  are  designated  as  slave  nodes  and  the  standard  elements  form  the  master  surface.  As  was 
the  case  for  the  attached  interface,  the  standard  nodes  (ma,  m3,  ns)  do  not  directly  affect  the  strain 
rates  in  SPH  node  i,  and  do  not  directly  affect  the  forces  on  node  i.  The  maximum  allowable 
overlap  of  an  SPH  node  and  slave  node  i,  with  master  surface  m2  -  m3,  is  shown  as  5o.  When  the 
SPH  nodes  are  initially  defined  as  SPH  nodes,  then  5o  =  0.  When  the  SPH  nodes  are  converted 
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Figure  33.  SPH  Node  Attachment  to  a  Standard  Grid 


from  standard  elements,  6o  can  be  a  small  fraction  of  the  SPH  node  diameter.  The 
approximations  for  the  sliding  interface  are  similar  to  those  of  the  attached  interface,  but  more 
research  is  required  to  evaluate  and  improve  the  existing  algorithm.  The  specific  equations  are 
presented  in  subsection  3.4.2. 

When  slave  node  i  overlaps  (5  >  5o)  master  segment  m2  -  m3,  then  the  three  normal  velocities  of 
nodes  i,  m2,  and  m3  are  adjusted  to  (1)  conserve  linear  momentum,  (2)  conserve  angular 
momentum  and  (3)  provide  a  normal  velocity  match  of  node  i  on  master  segment  m2  -  m3.  The 
positions  of  the  nodes  are  adjusted  to  be  consistent  with  the  velocity  changes.  The  specific 
equations  are  presented  in  subsection  3.4.2. 
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T19616 

Figure  34.  SPH  Node  Sliding  on  a  Standard  Grid 

It  is  also  possible  to  generate  SPH  nodes  from  a  standard  finite  element  grid.  Figure  8  shows  an 
illustration  of  a  standard  erosion  algorithm,  as  described  in  subsection  3.4.3,  and  Figure  35  shows 
the  SPH  node  generation  algorithm.  The  standard  algorithm  removes  (erodes)  elements  from  the 
sliding  interface  when  the  elements  are  strained  to  an  equivalent  plastic  strain  of  Ep  »  li ,  and  it 
fiien  redefines  the  master  surface  and  slave  nodes.  Although  the  mass  is  retained  at  the  nodes, 
the  volumes  of  the  eroded  elements  are  discarded,  and  this  introduces  inaccuracies  into  the 
algorithm. 

The  SPH  generation  algorithm  is  similar,  except  it  replaces  the  highly  strained  elements  on  the 
master  surface  with  equivalent  SPH  nodes.  Because  the  SPH  nodes  replace  both  the  mass  and 
volume  of  the  eroded  elements,  this  is  a  more  accurate  algorithm.  Here  the  elements  are 
converted  to  SPH  nodes  at  an  equivalent  plastic  strain  of  ~  0.5 . 

When  a  standard  triangular  element  is  replaced  by  a  circular  SPH  node,  the  circle  will  extend 
beyond  the  replaced  triangular  element  by  a  distance  5o,  as  shown  in  Figure  35.  For  the 
subsequent  sliding  interface  computations,  an  effective  crossover  distance  is  defined  as 
5efr  =  5  -  5o,  where  6  is  the  current  crossover  distance  and  5o  is  the  initial  crossover  distance. 
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Figure  35-  Description  of  the  SPH  Node  Generation  Aigorithm 

This  effective  distance  is  then  used  to  determine  contact  (5eff  ^  0),  and  to  adjust  the  velocities  and 
displacements  of  the  master  and  slave  nodes  as  discussed  previously. 

The  current  SPH  generation  algorithm  has  a  significant  approximation  in  that  the  generated  SPH 
nodes  are  allowed  to  slide  along  the  standard  elements,  instead  of  being  attached  to  the  elements. 
Future  effort  is  required  to  develop  an  interface  algorithm  that  will  allow  the  SPH  node  to  be 
attached  to  the  standard  finite  element  grid  after  it  has  been  generated. 


D  tig616.doc 


138 


SECTION  6 
MATERIAL  MODELS 


This  section  describes  the  material  models  that  are  available.  The  basic  material  types  are  as 
follows: 

•  Solid  (metal)  materials 

•  Explosive  materials 

•  Crushable/concrete  materials 

•  Liquid  materials 

•  Brittle  (ceramic)  materials 

•  Reactive  explosive  materials 

•  RDG  solid  (metal)  materials. 

The  solid  (metal)  material  models  and  algorithms  are  presented  in  detail  and  the  imique  portions 
of  the  others  are  also  presented.  In  some  instances,  there  are  many  models  available  for  a 
material  type,  and  in  other  cases  there  is  only  a  single  model. 

6.1  SOLID  MATERIALS 

The  solid  material  models  are  generally  used  for  metals,  although  some  of  the  models  can  be 
used  for  nonmetallic  materials.  The  stress  determination  algorithms  for  ID,  2D,  and  3D  are 
provided  in  subsections  3.1.3  (ID),  3.4.4  (2D),  and  3.8.4  (3D). 

The  3D  stresses,  as  presented  previously  in  subsection  3.8.4,  are  composed  of  three  normal 
stresses  (a„,ay,a2)  and  three  shear  stresses  ('t,y,x,j2,Ty2) .  The  three  normal  stresses  are 

expressed  as 


-(P  +  Q) 

(1) 

-(P  +  Q) 

(2) 

-(P+Q) 

(3) 
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where  Sx,  Sy,  and  Sz  are  the  normal  deviator  stresses,  P  is  the  hydrostatic  pressure,  and  Q  is  the 
artificial  viscosity. 

Trial  values  of  the  deviator  stresses  and  shear  stresses  at  time  =  t  +  At  are 


s^T""*  =s‘ +2Ge,At  +  As,  (4) 

=s;+2GeyAt  +  ASy  (5) 

s*^^' =s;+2Ge,At  +  As,  (6) 

T;^^=T*y+GYxyAt  +  Ax,y  (7) 

O  ='>^!cz  +  GYxzAt  +  ATx,  (8) 


In  Equation  4  the  first  term  (s^ )  is  the  normal  stress  at  the  previous  time  and  the  second  term 
(2GexAt)  is  the  incremental  stress  due  to  the  incremental  strain  (e^At)  during  that  time 
increment,  where  G  is  the  elastic  shear  modulus  and  At  is  the  integration  time  increment.  The 
third  term  (As^  )  is  due  to  shear  stresses  fi-om  the  previous  time  increment,  which  now  act  as 

normal  stresses  due  to  the  new  orientation  of  the  element  caused  by  an  incremental  rotation 
(©y  At,  ©xAt)  during  the  time  increment.  The  remaining  normal  stresses  and  shear  stresses  have 

a  similar  form. 

The  correction  terms  for  element  rotations  are  given  in  Equations  37—42  of  subsection  3.8.4. 

Equations  4—9  assume  an  elastic  response  of  the  material.  If  the  strength  of  the  material  is 
exceeded,  then  plastic  flow  (or  fracture)  will  occur.  The  Von  Mises  yield  criterion  is  used  to 
determine  an  equivalent  stress,  a ,  that  can  be  compared  to  the  uniaxial  tensile  (or  compressive) 
strength  of  the  material.  The  general  form  of  the  equivalent  stress  is 

O  =  -  a,  f  +  (o.  -  O.f  +  (o,  -  v' )]  ('“) 
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Using  deviator  stresses  (instead  of  total  stresses).  Equation  10  can  be  rewritten  as 


<7  =  -^|(sx  +  sj  +  s^)  +  +  t;,)  (11) 

If  o  is  not  greater  than  the  equivalent  tensile  strength  of  the  material,  o,  the  final  deviator  and 
shear  stresses  are  as  given  in  Equations  4-9.  If  a  is  greater  than  a,  then  the  stresses  in 
Equations  4-9  are  multiplied  by  the  factor  (g/ a) .  When  the  reduced  deviator  and  shear  stresses 
are  put  into  Equation  1 1,  the  result  is  always  g  =  g  .  This  is  known  as  the  radial  return 
algorithm. 

In  the  following  subsections,  models  are  presented  for  the  strength  of  the  material  (g),  the 
pressure  (P)  and  artificial  viscosity  (Q).  Additional  subsections  address  a  2D  orthotropic 
algorithm,  the  internal  energy  algorithm,  fracture  models,  and  a  fragmentation  model. 

6.1.1  Johnson-Cook  Strength  Model 

The  Johnson-Cook  strength  model  (Reference  37)  is  expressed  as 

G  =  [c,+C28^][n-C3lne*][l-T*"']-i-C4P  (1) 

Where  £  is  the  equivalent  plastic  strain,  e*  =  e  /  6^  is  the  dimensionless  strain  rate  for 

=  1.0s"’ ,  T*  is  the  homologous  temperature,  and  P  is  the  hydrostatic  pressure.  This  model  is 
valid  only  for  0  <  T*  <  1.0.  The  material  constants  are  Ci,  C2,  N,  C3,  M,  and  C4.  Ci  and  C2  have 
imits  of  stress  and  the  others  are  dimensionless.  The  original  model  did  not  include  the  pressure 
term  (C4P). 

Although  this  model  is  empirical,  it  is  flexible,  robust,  and  contains  the  effects  of  the  important 
parameters.  The  strength  goes  to  zero  as  the  temperature  approaches  melting  (T*  =  1.0).  It  is 
also  relatively  easy  to  obtain  constants  for  this  model  (References  37  and  38). 

A  constant  flow  stress  (g  =  Ci)  can  be  obtained  by  setting  C2  =  C3  =  C4  =  T**^  =  0,  a  linear  strain 
hardening  model  (g  =  Ci  +  C2  *  s)  can  be  obtained  by  setting  C3  =  C4  =  T**^  =  0,  and  a  pressure 
hardening  model  (g  =  Ci  +  C4P)  can  be  obtained  by  setting  C2  =  C3  =  T*’^  =  0. 
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6.1.2  Modified  Johnson-Cook  Strength  Model 

The  modified  Johnson-Cook  strength  model  (Reference  38)  is  expressed  as 

a  =  [c,+C2e''Ie*^’][l-T*"']+C4P  (1) 

This  is  identical  to  the  Johnson-Cook  model  in  subsection  6.1.1,  except  that  the  strain  rate  effect 
is  different.  This  model  provides  an  enhanced  strain  rate  effect  at  high  strain  rates. 

6.1.3  Zerilli'Armstrong  Strength  Model  for  FCC  Metals 

The  Zerilli- Armstrong  strength  model  (Reference  39)  for  FCC  metals  is  expressed  as 

a  =  Co  +  C2e''exp(-C3T+C4Tlne)  (1) 

Where  e  is  the  equivalent  plastic  strain,  T  is  the  absolute  temperature,  and  8  is  the  equivalent 

strain  rate.  The  material  constants  are  Co,  C2,  N,  C3,  and  C4.  Co  and  C2  have  units  of  stress,  and 

C3  and  C4  have  xmits  of  (temperature)"'.  The  grain  size  is  not  represented  as  a  variable  (as  it  is  in 

Reference  39)  but  is  included  in  Co- 

6.1.4  Zerilli-Armstrong  Strength  Model  for  BCC  Metals 

The  Zerilli-Armstrong  model  (Reference  39)  for  BCC  metals  is  expressed  as 

a  =  Co  +  C,  exp(-  C3T + C4T  Ine  *)-H  (1) 

This  is  similar  to  the  Zerilli-Armstrong  FCC  model  in  subsection  6.1.3.  Co,  Ci,  and  C5  have 

units  of  stress,  and  C3  and  C4  have  units  of  (temperature)"'. 

6.1.5  Bodner-Partom  Strength  Model 

The  Bodner-Partom  model  (References  40, 41,  and  42)  is  expressed  as 
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a  =  Z 


^  2N  ^ 

,N  +  lJ 


r 

•  In 

V 


(1) 


where  is  the  equivalent  plastic  strain  rate  and  Do  is  an  input  constant  representing  the 

maximum  allowable  plastic  strain  rate.  The  other  input  constants  are  Zo,  Zi,  No,  Ni,  Mo,  Mi,  and 
a. 


For  constants  Mi  =  a  =  0 

Z  =  Z,-(z,-Z„)exp(-M.W,) 


(2) 


where  Wp  is  the  plastic  work  per  unit  volume. 

For  constants  Mi  >  0  and  a  >  0  an  expanded  form  is  used. 

Z  =  Z,  -(z,  -Zo)exp(-MoWp)exp(-(Mo  +M,  -M)/a)  (3) 


where 

M  =  Mo  +  M,exp(-aWp) 


(4) 


ForNi  =Tzero 
N  =  No 

and  forNl  >  Tzero 

f N  _T  ^ 

N  =  No+  ' 

V  ^  ^zero  J 

where  T  is  the  absolute  temperature  and  Tzero  is  the  absolute  zero  temperature. 


(5) 


(6) 
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This  constitutive  model  requires  an  iterative  solution  because  the  strength,  a,  and  the  plastic 
strain  rate,  £p ,  are  dependent  upon  one  another.  This  results  in  increased  CPU  time.  Also,  this 

model  is  more  complex  than  the  Johnson-Cook  and  Zerilli- Armstrong  models,  and  References 
40, 41,  and  42  should  be  consulted  for  more  information. 


6.1.6  MTS  Strength  Model 


The  MTS  model  is  presented  in  References  43  and  44.  The  following  description  is  an  edited 
version  provided  by  P.J.  Maudlin  of  Los  Alamos  National  Laboratory.  The  material  strength  is 


+  s*,i 


(1) 


which  contains  the  constants  (which  represents  disloeation  interactions  with  long-range 
barriers),  d;  (which  represents  dislocation  interactions  with  interstitial  atoms,  and  (which 
represents  dislocation  interactions  with  solute  atoms). 

The  first  product  in  the  equation  for  a  contains  a  micro-structure  evolution  variable,  i.e.,  a , 

called  the  mechanical  threshold  stress,  that  is  multiplied  by  a  constant-structure  deformation 
variable  Sth;  Sth  is  a  function  of  absolute  temperature  T  and  plastic  strain-rate  Ep .  The  evolution 

equation  for  a  is  a  differential  hardening  law  representing  dislocation-dislocation  interaction: 


=  0. 


tanh 


'  a  ' 

a— 


1  — 


tanh(a) 


(2) 


where  is  the  value  of  d  at  large  plastic  strain  and  a  is  a  material  constant. 


In  the  equation  for  ^ — ,  ©q  represents  hardening  due  to  dislocation  generation  and  the  stress 
d£p 

ratio  represents  softening  due  to  dislocation  recovery.  The  threshold  stress  at  zero  strain¬ 
hardening  dj  is  called  the  saturation  threshold  stress.  Relationships  for  Og ,  are: 


0.  =an-i-a 


,ln(ep)+a2.J^ 


(3) 
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which  contains  the  dislocation  generation  material  constants  ai,  az,  and  as. 
Also 


\kT/Gbl 


v^soy 


(4) 


where  6^  is  the  saturation  stress  at  zero  degrees  kelvin,  e^o  is  a  reference  strain  rate,  b  is  the 
magnitude  of  Btirgers  vector  (inter-atomic  slip  distance),  A  is  another  material  constant,  and  K  is 
Boltzmanns  constant. 


The  shear  modulus  in  these  equations  is  assumed  to  be  a  function  of  the  temperahare  and  is 
expressed  as 


G  =  bo-b, /(exp(b2T)-l)  (5) 

where  bo  is  the  shear  modulus  at  zero  degrees  Kelvin,  and  bi  and  hz  are  constants. 

For  thermal-activation  controlled  deformation  Sth  is  evaluated  via  an  Arrhenius  rate  equation: 

I 
p 

(6) 


S.K  = 


1- 


kTln(eo/ep) 


Gb^g„ 


where  e,,  is  a  reference  strain  rate,  go  is  the  normalized  activation  energy  for  a  dislocation 
interaction,  and  p  and  q  are  additional  constants. 

Expressions  for  Sth,i  and  Sth,s  are  identical  to  the  equation  for  Sth  in  form  but  use  the  constants  e,,  j , 
go,i, pi,  qi  (for  Sth.i)  and  So,,  go,s, Ps,  qs (for  Sth,s). 

The  MTS  model  is  much  more  complex  than  the  Johnson-Cook  and  Zerilli-Armstrong  models, 
and  References  43  and  44  should  be  consulted  for  more  information. 
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6.1.7  Steinberg-Guinan-Lund  Strength  Model 

The  Steinberg-Guinan-Lund  models  are  described  in  References  45  and  46.  A  rate  independent 
model  is  given  in  Reference  45  and  a  rate  dependent  model  is  provided  in  Reference  46.  In  this 
subsection  the  material  strength  is  represented  by  Y  instead  of  a. 


The  strain-rate  dependent  form  of  the  SGL  model  defines  the  yield  stress  as 


Y 


Y,(£„T)+Y,f(e,)] 


G(P,T) 

Go 


(1) 


where  the  athermal  and  thermally  activated  components  Ya  f  (Sp)  and  Yt  are  defined  by 


Y»f(e,)  =  Y4l+p(e, +£,)}" 


(2) 


and 


£ 


p 


f 

\ 


2Uk 

T 


1 


(3) 


In  these  equations,  P  and  T  are  the  pressure  and  temperature,  Ep  and  are  the  equivalent  plastic 

strain  and  equivalent  plastic  strain  rate,  Ya  is  the  yield  strength  at  the  Hugoniot  elastic  limit 
(HEL)  and  f(  Ep )  is  the  work-hardened  function  with  (P,  6i,  and  n}  as  fitting  parameters.  Go  is 

the  initial  shear  modulus,  Yp  is  the  Peierls  stress,  and  2Uk  is  the  energy  necessary  to  form  a  pair 
of  kinks  in  a  dislocation  segment.  The  quantities  Ci  and  C2  are  defined  in  terms  of  various 
dislocation  mechanics  parameters  and  are  specific  to  the  material  being  modeled.  The  shear 
modulus  G  is  defined  by 


G(P,T)  =  Go 


1  ^1/3  Xooro ) 


(4) 


with  A  and  B  treated  as  material  constants  and  q(p/po)  denoting  the  compression.  Melting  is 
modeled  using  a  modified  Lindemann  law,  where  the  melt  temperature  Tm  is  defined  as 
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T™=T™oexp 


f  01 

2a 

1  — 

1  ^J] 

11' 


2(Yo-a-l/3) 


(5) 


Tmo  is  the  melt  temperature  at  constant  volume,  yo  is  the  initial  Gruneisen  coefficient,  and  a  is  a 
material  constant.  For  any  compression  r\  in  which  the  temperature  T  exceeds  Tm,  melting  is 
considered  to  have  occurred,  resulting  in  the  loss  of  yield  strength  and  shear  strength.  In 
addition,  there  are  two  limits  imposed. 

Y.f(£,)sY;„  (6) 

And 

Yt<Yp  (7) 

where  is  the  work-hardening  maximum  in  the  rate-dependent  version  of  the  model.  The 
rate-dependent  form  of  the  SGL  model  is  treated  as  a  special  case  of  the  more  general  form  for 
yield  stress  with  Yj  set  to  zero.  In  particular,  the  rate-independent  model  assumes  the  form 

Y  =  Yof(ep)G(P,T)/Go  (8) 

where  f(ep)  and  G(P,T)  are  defined  as  before.  Here,  however,  the  following  limit  applies: 

V(ep)sY^  (9) 

These  models  are  more  complex  than  the  Johnson-Cook  and  Zerilli-Armstrong  models,  and 
References  45  and  46  should  be  consulted  for  more  information. 

6.1.8  2D  Orthotropic  Model 

This  is  a  simplified  anisotropic  model  that  is  anisotropic  only  in  the  plastic  response.  The  elastic 
response  remains  isotropic.  It  essentially  expands  and  contracts  the  yield  surface  in  the  three 
primary  axes  and  allows  for  initial  orientation  of  the  material  axes  with  the  system  coordinate 
axes.  It  is  described  in  detailed  in  Reference  47.  The  description  that  follows  is  an  edited 
version  provided  by  P.J.  Maudlin  of  Los  Alamos  National  Laboratory. 
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The  1948  Hill  yield  function  (Reference  48)  written  in  terms  of  total  stress  (ay)  is 

+H(aj]— 022) 

+ 2Naj2  +  2Maj3  +  2La^  }  -  =  0 

where  the  quantity  a  is  a  flow  stress  that  is  assumed  to  be  a  function  of  strain,  strain  rate,  and 
temperature  invariants.  Expressing  Hill’s  Quadratic  in  terms  of  deviatoric  stress  (sy)  gives 

f  =  (F"*"  H)S22  +(F+  G)S33  —  2HS]iS22  “2Gs,|S33  —  2FS22S23 

+  2Ns^2  +  2Ms^3  +  2LS23  }  -  s^  =  0 


The  physical  interpretation  and  measurement  of  the  parameters  in  Equation  2  are  discussed  by 
Hill  in  some  detail.  This  function  is  conveniently  analytic  and  can  be  used  for  three-dimensional 
problems.  In  terms  of  independent  stress  components  it  represents  a  five-dimensional  stress 
space. 

For  a  two-dimensional  problem,  recalling  that  the  trace  of  the  deviatoric  stress  tensor  is  zero. 
Hill’s  yield  function  reduces  to 

f  =  I  {(4G  -H  H  +  F)s^,  +  (4F  4-  H  +  G)s^2  +  (4F + 4G  -  2H)s„S22  -h  2N  s  sf2 }  -  a'  =  0  (3) 

The  yield  surface  actually  implemented  in  EPIC  is  a  u-plane  transformation  (useful  for  yield 
surface  data  fitting)  of  Equation  3  defined  in  terms  of  three  new  variables  (Sx,  Sy,  Sz)  which  are 
linear  combinations  of  the  deviatoric  stress  components. 


> 

Sx 

r-V3 

2 

-V3 

0 

3 

0 

Sy 

2 

0 

S22 

0 

1 

0 

-V3 

J 

These  transformed  stresses  represent  an  orthogonal  reference  fi:ame  whose  Sx  -  Sy  plane  is 
coincident  with  the  7t-plane  in  the  principle  axes  space  for  deviatoric  stress.  Substituting 
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Equation  4  into  Equation  3,  assuming  that  the  rotational  term  SxSy  is  insignificant,  and  collecting 
terms,  gives  the  form  that  is  used  in  EPIC  for  two-dimensional  problems. 


f  =  As^-t-BsJ-i-Cs^-a^=0 


(5a) 


Or 


f  =  Af-fs„-V3: 


»22 


Y 

f3  ^ 

+  B 

— s,, 

J 

U  "j 

-hC(-V3s,2)'-a^=0 


(5b) 


where  the  shape  factor  constants  A,  B,  and  C  can  be  constructed  from  the  Hill  constants  F,  G,  H, 
and  N.  For  the  special  case  where  the  shape  factors  A  =  B  =  C  =  1,  Equation  5  simplifies  to  the 
von  Mises  yield  fimction. 

The  yield  fimction  given  by  Equation  5  may  contain  an  implied  normalization  in  the  quantity  cr. 
As  it  appears  in  Equation  5,  ct  represents  a  directionally  averaged  flow  stress.  If  the  flow  stress  is 
constructed  solely  from  \miaxial  stress  data  in  a  given  material  direction,  such  as  “1”  for 
example,  then  the  flow  stress  function  Gus  (in  equivalent  stress  imits)  is  conceptually  different 
fi’om  the  quantity  appearing  in  Equation  5  and  needs  to  be  renormalized,  i.e., 

a  =  Ous/M  (6) 

in  order  to  recover  the  imiaxial  stress  result: 

Su=±fous  (7) 


Substitution  of  Equations  6  and  7  into  the  yield  function  given  by  Equation  5  and  solving  for  the 
normalization  constant  M  under  the  assumption  of  uniaxial  stress  in  the  “1”  direction  gives 


The  important  point  here  is  that  if  a  flow  stress  fimction  like  Johnson-Cook  is  used  for  a  in 
Equation  5,  where  the  model  coefficients  have  been  characterized  fi’om  data  measured  in  a 


certain  direction,  then  renormalization  of  the  flow  stress  model  by  some  M  is  necessary. 
Obviously  as  the  yield  function  tends  to  isotropy,  Equation  8  indicates  that  M  goes  to  unity. 

The  yield  function  given  by  Equation  5  is  cast  in  a  reference  frame  that  rotates  and  translates  with 
an  element  of  material  relative  to  the  laboratory  frame.  From  the  point  of  view  of  an  observer  in 
this  material  frame,  the  material  element  experiences  only  deformation  without  rigid  body 
motion  (pure  stretch).  For  two-dimensional  problems  the  rigid  body  rotations  are  planar 
occurring  in  the  x-z  plane,  and  are  described  by  the  angle  0.  An  initial  value  for  this  angle,  i.e., 

0  (t  =  0),  can  be  specified  that  measures  the  initial  orientation  difference  between  the  z  axis  in  the 
laboratory  frame  and  the  z  axis  in  the  material  frame.  Thus,  the  material  constants  for  Equation  5 
(i.e..  A,  B,  C)  can  be  measured  in  a  reference  frame  initially  different  from  the  reference  frame 
that  EPIC  uses  for  some  application  calculation. 

6.1.9  Internal  Energy 

The  shear  and  deviator  stresses,  together  with  the  shear  and  deviator  strain  rates,  generate  internal 
energy  in  the  material.  This  energy  leads  to  increased  temperatures  (that  affect  the  strength 
models)  and  it  is  also  used  in  the  Equation  of  State  model  for  the  pressures.  This  incremental 
internal  energy  (for  a  cycle  of  integration)  is 

AE,  =  (s,e,  +  s^h^  +^xztxz  +  VYyz)  (v/ vj  At  (1) 

The  bars  on  the  deviator  and  shear  stresses,  and  the  volume,  represent  averages  of  tiiese  values  at 
times  t  and  t  +  At.  The  factor  V  /  V„  converts  the  energy  to  internal  energy  per  initial  volume. 


6.1.10  Artificial  Viscosity 

The  artificial  viscosity  is  an  important  term  in  the  representations  for  the  normal  stresses,  as 
shown  in  Equations  1-3  in  subsection  6.1.  It  is  combined  with  the  normal  stresses  to  damp  out 
localized  oscillations  of  the  concentrated  masses.  It  tends  to  eliminate  spurious  oscillations 
which  would  otherwise  occiu"  for  wave  propagation  problems.  This  technique  was  originally 
proposed  by  Von  Neumann  and  Richtmyer  (Reference  33)  and  has  been  expanded  for  use  in 
various  eomputer  codes  (Reference  5).  It  is  expressed  in  terms  of  linear  and  quadratic 
components  and  is  applied  only  when  the  volumetric  strain  rate  (e^)  is  negative. 
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Q  =  CLpc,h|eJ  +  Cgph^  (e fore,  <0 


(1) 


Q  =  0  for  e,  >  0 

where  c*  is  the  soxmd  velocity  of  the  material,  p  is  the  density,  and  h  is  a  characteristic  dimension 
(such  as  minimum  altitude)  of  the  element.  Cl  and  Cq  are  the  linear  and  quadratic  dimensionless 
coefficients. 

6.1 .1 1  Mie-Gruneisen  Equation  of  State 

The  Mie-Gruneisen  Equation  of  State  (Reference  49)  is  used  to  compute  the  pressure  term  in 
Equations  1-3  of  subsection  6.1.  The  pressme  is  a  function  of  the  volumetric  strain  and  the 
internal  energy. 

p  =  p,+rE.(i+n)  (1) 

where 

P,  =  (k,p  -H  ^  j  (2) 


Substituting  Equation  2  into  Equation  1  gives  the  complete  expression  for  the  pressure. 


P  =  (K,^  +  +  K,n=)  [^1  +  rE,(H- ^) 


(3) 


where  p  =  V(/V-l,  Ki,  K2,  and  K3  are  material-dependent  constants  and  F  is  the  Gruneisen 
coefficient.  Vo  and  V  are  the  initial  and  current  element  volumes. 


Since  the  pressure  can  be  significantly  affected  by  the  internal  energy  Es,  it  is  desirable  to  solve 
the  pressure  and  energy  equations  simultaneously.  This  gives 


Et+A.  ^ 


E‘-  .5[(P-HQ)‘+Q*^'“-l-Pr'^ 


£^At  H-  AEj 


1  -I-  .5^(l-l-^)e,At 


(4) 
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where  is  the  volumetric  strain  rate  and  AEa  is  the  internal  energy  generated  by  the  deviator 
and  shear  stresses  during  the  previous  cycle.  Es  is  the  total  internal  energy  per  initial  imit 
volume. 

Substituting  Equation  4  into  Equation  3  gives  the  pressure  at  time  =  t  +  At  for  the  internal  energy 
at  time  =  t  +  At. 

It  is  interesting  that  Equation  3  reduces  to  P  =  Kip  for  small  strains  and  spherical  stress  (no 
shear).  The  higher  order  terms  and  Ksp^)  go  to  zero  and  the  internal  energy  is  K,pV2.  This 
shows  that  Kj  is  consistent  with  an  elastic  buUc  modulus. 


Another  form  of  the  Mie-Gruneisen  Equation  of  State  is 


P  PqCsTI 

■(1-snr 


+rE. 


(5) 


where  Po  is  the  initial  density,  q  =  1  -  Po/p,  E,  is  the  internal  energy  per  initial  vmit  volume,  c,  is 
the  bulk  sovmd  velocity,  T  is  the  Gruneisen  coefficient,  and  s  is  the  slope  of  the  U*  -  Up 
relationship.  -  Up  are  the  soxmd  velocity  and  particle  velocity,  respectively.  Again  the 
pressure-energy  equations  are  solved  simultaneously  as  they  are  in  Equations  1—4. 

6.1.12  Johnson-Cook  Fracture  Model 

The  Johnson-Cook  firacture  model  (Reference  50)  is  based  on  accumulated  damage.  When  the 
damage  approaches  imity  (D  =  1 .0),  the  element  has  failed.  After  the  element  has  failed  the 
material  essentially  behaves  as  a  liquid  because  it  has  no  strength  (no  shear  and  deviator  stresses) 
and  it  cannot  develop  hydrostatic  tension.  A  failed  element  can  produce  only  hydrostatic 
compression.  It  is  also  possible  to  soften  (weaken)  the  strength  as  the  damage  increases,  such 
that  the  failure  is  gradual  rather  than  instantaneous. 


The  damage  to  an  element  is  defined  as 


where  Asp  is  the  increment  of  equivalent  plastic  strain  which  occurs  during  an  integration  cycle, 
Ep  is  the  equivalent  strain  to  fracture,  under  the  current  conditions  of  strain  rate,  temperature, 

pressure,  and  equivalent  stress.  Fracture  is  then  allowed  to  occur  when  D  =  1 .0. 

To  illustrate,  assume  the  fracture  strain  in  Figure  36  is  a  function  of  only  the  pressure-stress  ratio, 
a*,  which  will  be  defined  later.  This  assumed  relationship  indicates  the  material  can  be  more 
severely  strained  when  it  is  under  hydrostatic  compression.  The  example  shows  how  damage  is 
accumulated  under  compression.  When  the  pressure  is  released  to  a*  =  0,  D  =  0.75  even  though 
the  plastic  strain,  gp,  is  greater  than  at  a*  =  0.  The  importance  of  path-dependency  is  evident 

since  a  model  based  only  on  current  conditions  would  incorrectly  predict  fracture  for  these 
conditions. 


Assumed  Relationship  Example 


Figure  36.  Description  of  the  Johnson-Cook  Fracture  Modei 

The  general  expression  for  the  strain  at  fracture  is  given  by 

ej  =[d,  +D2  exp'^^'^'ll+D^  InE^ll-nD^T*]  (2) 

for  constant  values  of  the  variables  (<t*,  e  * ,  T*)  and  for  ct*  <  1 .5.  The  dimensionless  pressure- 
stress  ratio  is  defined  as  a*  =  cy„  /  a  where  Om  is  the  average  of  the  three  normal  stresses  and  a 
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is  the  von  Mises  equivalent  stress.  The  dimensionless  strain  rate,  e  *,  and  homologous 
temperature,  T*,  are  identical  to  those  used  in  the  strength  models  presented  in  the  previous 
subsections. 


The  five  dimensionless  constants  are  Di  ...  D5.  The  expression  in  the  first  set  of  brackets 
follows  the  form  presented  by  Hancock  and  Mackenzie  (Reference  51).  It  essentially  says  that 
the  strain  to  fi-acture  decreases  as  the  hydrostatic  tension,  am,  increases.  The  expression  in  the 
second  set  of  brackets  represents  the  effect  of  strain  rate,  and  that  in  the  third  set  of  brackets 
represents  the  effect  of  temperature.  For  high  values  of  hydrostatic  tension  (a*  >  1.5)  a  different 
relationship  is  used. 

It  is  clear  from  Figure  36  that  less  damage  is  accumulated  when  the  material  is  in  compression 
(a*  <  0).  At  the  other  extreme,  however,  when  there  is  significant  hydrostatic  tension,  the  strain 
to  fi-acture  drops  rapidly.  Figure  37  shows  the  relationship  which  is  used  for  large  values  of  the 
pressure-stress  ratio  (a*  >  1.5).  The  fiacture  strain  varies  in  a  linear  manner,  fiom  a*  =  1.5  to 
a[pa„  at  .  The  model  constants  are  the  spall  stress,  aspaii,  and  the  minimum  fiacture  strain, 

.  The  dimensionless  is  computed  fiom  Cspaii  and  the  current  value  of  the  von  Mises 
flow  stress,  a . 


6.1.13  Modified  Johnson-Cook  Fracture  Model 


The  modified  Johnson-Cook  fiacture  model  is  similar  to  the  Johnson-Cook  fiacture  model 
(presented  in  subsection  6.1.12)  except  the  Tuler-Butcher  time-dependent  spall  model  (Reference 
52)  is  used  for  the  high  tensile  regions.  For  mean  tensile  pressures  less  than  a™,  the  damage  and 
fracture  strain  are  determined  fiom  Equations  1  and  2  of  subsection  6.1.12.  For  mean  tensile 
stresses  a^  >  the  damage  is 


D  = 


IK-l^At 


K' 


(1) 


where  ,  At  is  the  integration  time  increment,  and  X  and  K*  are  material  constants. 
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Description  of  Fracture  Strains  at  Large  Tensiie  Pressure-Stress  Ratios 
for  the  Johnson-Cook  Fracture  Model 


6.1.14  Fragmentation  Model 

The  fragmentation  model  is  presented  in  Reference  18.  The  theoretical  approach  for  this  work  is 
an  extension  of  work  developed  by  Grady  (References  53, 54,  and  55)  on  the  fragmentation  of 
stretching  jets  and  expanding  sheets.  Fragment  sizes  for  one-dimensional  stretching  jets  and 
two-dimensional  expanding  sheets  were  developed  using  energy  principles.  These  theories  are 
based  on  a  ductile  mode  of  failure.  Fragment  sizes  are  determined  by  equating  the  internal 
kinetic  energy  of  the  flowing  fragment  material  with  the  plastic  strain  energy  used  during  the 
fragmentation  process.  The  following  discussion  takes  Grady’s  theory  and  extends  it  to  three- 
dimensional  ductile  fragmentation,  with  considerations  being  taken  to  ensure  computational 
compatibility. 

Consider,  prior  to  fragmentation,  a  rapidly  flowing  mass  of  material  from  which  a  cubical 
fragment  of  size  b  will  be  formed  as  shown  in  Figure  38.  To  simplify  the  development,  assume 


DJ19616.doc 


155 


that  the  x,  y,  and  z  axes  correspond  to  those  of  the  principal  strain  rates.  The  theoretical 
approach  is  to  equate  the  kinetic  energy  of  a  fragment’s  flowing  material  about  its  center  of  mass 
(local  kinetic  energy),  witih  the  plastic  work  occurring  during  the  fragmentation  process.  If  it  is 
assumed  that  all  of  the  internal  kinetic  energy  is  consumed  during  fragmentation,  an  equality  can 
be  developed  and  solved  directly  for  the  firagment  size. 


Figure  38.  Fragment  Characteristics 

The  local  kinetic  energy  for  the  flowing  fragment  material  can  be  determined  explicitly  by 
Slimming  the  local  kinetic  energy  in  each  direction: 

KE  =  -^J  v(x)^dm+^J  v(y)  dm+^J  v(z)  dm  (1) 

The  local  kinetic  energy  in  the  x-direction,  about  the  fragment’s  center  of  mass  can  be  expressed 


(KE),= 


(2) 
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where  the  incremental  mass  is 


dm  =  p  dx 

and  where  p  is  the  density  of  the  material. 

For  a  linear  velocity  gradient  across  the  fragment,  with  local  face  velocities  of  ±  Vx  (about  the 
CG  velocity),  the  x  velocity  distribution  is 


f  2x 

v(x)  =  vJ— -1 


forO<x<b 


The  local  kinetic  energy  in  the  x-direction  is  then 


fb  1  -  2(2x  Y 

(KE),=(-bVv“[--lJdx 


b^p  , 

(KE), 


Since  v(x)  is  linear,  the  deviator  strain  rate,  e,^ ,  is  constant  within  the  fragment  and  can  be 
expressed  as  a  function  of  the  velocity 


The  local  kinetic  energy  in  terms  of  the  de  viator  strain  rate  and  fragment  size  is  then 


(KE),= 


b*pe^ 


The  same  approach  is  \ised  to  determine  the  local  kinetic  energy  in  the  y  and  z  directions,  which 
gives  a  total  local  kinetic  energy  of  the  fragment  of 
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The  deviator  strain  rates  can  be  expressed  in  terms  of  an  equivalent  strain  rate 


and  the  total  kinetic  energy  becomes 


(10) 


(11) 


The  plastic  work  generated  for  each  fragment  during  the  fragmentation  process  can  be  expressed 
in  the  following  form: 

PW  =  b=(s.e:+s,e;+s,e:)  (12) 

where  Sx,  Sy,  and  Sz  are  constant  deviator  stresses  and  ej ,  and  are  the  deviator  strains 

occurring  from  void  initiation  to  material  separation. 

Equation  12  can  be  expressed  in  terms  of  equivalent  stress  and  strain  to  obtain  the  following 
relationship  for  plastic  work: 

PW  =  b^aef  (13) 

Equation  13  represents  the  total  plastic  work  needed  for  fragmentation  to  occur,  where  is 
analogous  to  the  equivalent  strain  rate  of  Equation  10,  and  the  von  Mises  equivalent  stress  is 


Equating  Equations  1 1  and  13  and  solving  for  b,  the  fragment  size  is  expressed  as 
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The  dimensionless  constant,  A  =  ,  can  be  approximated  by  establishing  bounds  for 

Grady  (Reference  54)  estimated  that  the  strains  which  occur  from  void  initiation  to  material 
separation  are  between  0.01  and  0.10.  Substituting  0.01  <  Ef  <  0.10  gives  0.40  <  A  <  1.26. 


Equation  15  is  expressed  in  terms  of  equivalent  strain  rate  and  equivalent  stress  for 
computational  compatibility.  As  applied  munerically,  the  fragment  size  calculation  takes  the 
following  form: 


b  = 


(16) 


where  b  is  the  fragment  size  computed  for  the  current  element  condition  of  a,  p,  and  e  from 
Equation  15;  and  Afip  is  the  increment  of  plastic  strain  occurring  during  an  integration  cycle. 
The  total  plastic  strain  is  represented  by  Ep  =  ZAEp . 


This  calculation  is  performed  for  every  element  at  each  integration  cycle  until  the  element  fails. 
Here  the  failure  does  not  refer  to  the  fragmentation  presented  herein,  but  to  fracture  due  to 
damage  or  erosion.  Under  both  of  these  conditions,  the  equivalent  stress  is  set  to  a  =  0 . 

When  an  element  fails,  a  fragment  of  size  b  is  saved  in  an  element  array  for  post-processing. 
The  calculated  fragment  size  is  an  average  based  on  the  amoimt  of  element  strain  that  occurs  for 
each  cycle  as  shown  in  Equation  16. 

This  is  a  deviation  from  the  theoretical  approach  that  states  the  calculation  should  start  at  the 
initiation  of  fragmentation  and  continue  until  complete  separation.  Due  to  the  difficulty  in 
determining  when  the  initiation  of  fragmentation  begins,  the  approach  here  is  to  calculate  an 
average  firagment  size  based  on  the  element  history,  and  to  later  apply  a  factor  A  to  the  fragment 
size  during  post-processing.  Equation  16  is  the  only  addition  to  the  main  routine  in  EPIC  which 
is  required  to  calculate  the  fragment  size,  and  die  increased  computational  time  is  negligible. 

The  remainder  of  the  firagment  computations  are  performed  in  the  Postprocessor.  The  final 
firagment  size  associated  with  each  failed  element  is  obtained  by  applying  the  factor  A  to  the 
calculated  value  determined  from  Equation  16,  during  post-processing. 
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The  number  of  fragments  for  eaeh  failed  element  is  determined  by  dividing  the  element  volume 
by  the  fragment  volume.  Mass  is  conserved,  and  thus,  the  number  of  fragments  calculated  will 
generally  not  be  a  whole  number.  A  velocity  vector  is  defined  for  the  fragments  by  using  an 
average  of  the  nodal  velocities  associated  with  the  element.  All  the  fragments  that  are 
determined  from  a  specific  failed  element  will  have  the  same  velocity  vector.  For  each  failed 
element,  a  fragment  size,  velocity,  and  number  of  fragments  will  be  known.  The  fragment  mass, 
kinetic  energy,  and  momentum  are  determined  straightforwardly. 

The  EPIC  Postprocessor  can  plot,  in  bar  graph  form,  the  number  of  fragments  versus  fragment 
mass,  momentum,  kinetic  energy,  and  size.  The  advantage  of  presenting  the  data  in  this  form  is 
that  the  results  are  relatively  mesh  independent. 

6.2  EXPLOSIVE  MATERIALS 

The  materials  identified  as  Explosive  Materials  are  designated  to  detonate  at  a  predetermined 
time,  based  on  the  point  (or  points)  of  detonation  and  the  distance  from  the  initial  detonation 
point  (Reference  19).  There  are  also  Reactive  Explosive  Materials  (described  in  subsection  6.6) 
that  do  not  detonate  at  a  specified  time,  but  rather  are  detonated  when  the  proper  conditions 
(generally  the  pressure  history)  have  been  experienced. 

6.2.1  Gamma  Law  Model 

The  Gamma  law  explosives  model  determines  the  pressure  from 

p  =  f(y-i)e/v  (1) 

where  F  is  the  bum  fraction  (0  <  F  <  1 .0),  y  is  a  material  constant,  and  E  is  the  internal  energy  per 
initial  unit  volume.  The  relative  volume  is  V  =  V  /  V„ ,  where  V  and  V„  represent  the  current  and 

initial  element  volumes,  respectively. 

The  explosive  is  effectively  initiated  with  the  bum  fraction,  which  is  dependent  on  the  time  for 
the  detonation  wave  to  arrive  and  travel  through  the  element,  or  the  compressed  state  of  the 
element.  The  bum  fractions  for  these  two  conditions  are 
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(2a) 


F  = 


(t-tjD  +  b/2 
b 


F  = 


1-V 

i~ycj 


(2b) 


In  Equation  (2a),  t  is  the  current  time  and  ts  is  the  time  required  for  the  detonation  wave  to  reach 
the  center  of  the  element  when  traveling  at  the  detonation  velocity,  D.  A  reference  distance,  b,  is 
used  to  spread  the  wave  front  over  a  limited  number  of  elements.  For  ID  elements  b  =  Azo 

(where  Azo  is  the  initial  length  of  the  element),  for  2D  triangular  elements  b  =  3^/a7  (where  Ao 

1 

is  the  initial  area  of  the  element),  and  for  3D  tetrahedral  elements  b  =  2  VqS  (where  Vo  is  the 
initial  volume  of  the  element).  Equation  (2b)  gives  the  bum  fraction  in  terms  of  the  compressed 
state,  where  Vcj  =  y/(y  +  1)  is  the  Chapman- Jouquet  relative  volume.  This  allows  a  converging 
detonation  wave  to  travel  at  a  velocity  greater  than  D.  The  maximum  value  of  F  from  Equations 
(2a)  and  (2b)  is  selected,  if  it  is  within  the  limits,  0  <  F  <  1 .0.  If  F  is  negative  or  greater  than 
unity,  then  F  is  set  to  0  or  1.0,  respectively. 

Because  the  pressure,  P,  is  directly  proportional  to  the  internal  energy,  E,  the  pressme-energy 
equations  are  solved  simultaneously  in  a  manner  similar  to  that  used  for  the  Mie-Gruneisen 
Equation  of  State  in  subsection  6.1.1 1. 

It  is  possible  to  determine  y  as  a  function  of  the  other  material  parameters. 

Y  =  Vi  +  D'Po/2Eo  (3) 

where  D  is  the  detonation  velocity,  po  is  the  initial  density,  and  Eo  is  the  initial  internal  energy  per 
initial  xmit  volume. 

6.2.2  JWL  Model 

The  JWL  explosives  model  (Reference  56)  also  is  commonly  used.  Here  the  pressure  is 
expressed  as 

P  =  f[c,(i-C5/C2V)  exp  (-C2V)-i-C3(l-C5/C4V)  exp  (-C4V)-hC5E/ v]  (1) 
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where  F  is  the  bum  fraction  from  subsection  6.2.1,  E  is  the  internal  energy  per  initial  unit 
volume,  and  V  =  V  /  Vo  is  the  relative  volume.  C,,  C2,  Q,  C4,  and  C5  are  material  constants. 

Again  the  pressure-energy  equations  are  solved  simultaneously. 

6.3  CRUSHABLE/CONCRETE  MATERIALS 

There  are  two  models  available  for  crushable/concrete  materials.  Both  models  were  developed 
for  concrete,  but  they  also  can  be  used  for  other  cmshable  materials.  For  the  solid  (metal) 
materials  under  subsection  6.1,  it  is  possible  to  use  combinations  of  various  strength,  pressure, 
and  fracture  models.  For  the  concrete  models,  however,  the  strength,  pressure  and  damage  act 
together  to  form  a  single  model. 

6.3.1  Modified  Osborn  Model 

The  Osborn  concrete  model  (Reference  57)  is  composed  of  a  pressure-hardening  strength  model 
and  a  cmshable  pressure  model.  The  strength  is  expressed  as 

a  =  [c,(l-D)-HC4P][l-i-C3lne*]  (1) 

The  softening  due  to  damage  (1  —  D)  and  the  strain  rate  effect  [l  +  C3  Ine  *]  are  modifications 
made  to  the  original  Osborn  model.  It  is  recommended  that  the  damage  model  not  be  used,  and  it 
is  therefore  not  described  herein. 

In  Equation  1,  D  is  the  damage  (not  used  and  set  to  D  =  0),  P  is  the  pressure,  and  e*  =  e  /  e„  is 
the  dimensionless  strain  rate  for  gj,  =  1.0  s”’ .  The  three  constants  are  Ci,  C4,  and  C3,  where  Ci 
has  units  of  stress  and  the  other  (C4  and  C3)  are  dimensionless. 

The  model  for  the  pressure  is  shown  in  Figure  39.  The  constants  are  Pcmsh,  Pcmsh,  Ki,  K2,  K3, 
Kiock,  and  piock  as  shown  in  Figure  39.  The  maximum  hydrostatic  tension  is  Pmin-  For  increasing 
p  the  material  behaves  in  a  linear  elastic  manner  for  p  <  pcmsh-  The  transition  region  is  for 
Pcrush  <  P  <  Piock,  where  the  pressure  is  expressed  as 

P  =  Pcrush  +  KiP  +  Kjp'  -h  K3P:'  (2) 

where  p  =  It  -  Pcmsh  •  P  -  Piock  the  pressure  is  determined  from  Kiock  as  shown  in  Figure  39. 


7 


Figure  39.  Description  of  the  Pressure-Volume  Relationship  for  the  Osborn 
Concrete  Model  (Specific  Data  Shown  are  for  a  Typical  Concrete) 

For  decreasing  p  in  the  transition  region,  the  pressure  unloads  along  an  interpolated  bulk 
modulus  (between  elastic  and  locked). 

For  the  pressure  model,  Pcwsh,  Ki,  K2,  K3,  Kiock,  and  Pmin  have  imits  of  stress  (or  pressure),  and 
the  others  (pcrash  and  piock)  are  dimensionless. 
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6.3.2  Holmquist-Johnson-Cook  Model 

The  Holmquist-Johnson-Cook  concrete  model  (Reference  58)  is  summarized  in  Figure  40.  The 
strength  portion  of  the  model  is  shown  at  the  top  of  Figure  40.  The  normalized  equivalent  stress 
is  defined  as  a*  =  o  /  f ' ,  where  a  is  the  actual  equivalent  stress  and  f '  is  the  quasi-static 

xmiaxial  compressive  strength.  The  specific  expression  is 

a*  =  [A(l-D)  +  BP*^][l  +  Clne*]  (1) 


where  D  is  the  damage  (0  ^  D  ^  1.0),  P*  =  P  /  f'  is  the  normalized  pressure  (where  P  is  the  actual 
pressure),  and  e*  =  e  /  £„  is  the  dimensionless  strain  rate  (where  e  is  the  actual  strain  rate  and 
=  1.0s~'  is  the  reference  strain  rate).  The  normalized  maximum  tensile  hydrostatic  pressure  is 
T*  =  T  /  f',  where  T  is  the  maximum  tensile  hydrostatic  pressure  the  material  can  withstand. 

The  material  constants  are  A,  B,  N,  C  and  a  max,  where  A  is  the  ct  intercept  at  P*  =  0,  B  is  the 
normalized  pressure  hardening  coefficient,  N  is  the  pressure  hardening  exponent,  C  is  the  strain 
rate  coefficient,  and  is  the  normalized  maximum  strength  that  can  be  developed. 

The  damage  for  fracture  is  shown  in  the  lower  left  comer  of  Figure  40.  It  is  accumulated  from 
both  equivalent  plastic  strain  and  plastic  volumetric  strain,  and  is  expressed  as 


„Ae,+An„ 


(2) 


where  A£p  and  App  are  the  equivalent  plastic  strain  and  plastic  volumetric  strain,  respectively, 
during  a  cycle  of  integration;  and  £p  +  Pp  =  f(P)  is  the  plastic  strain  to  fracture  under  a  constant 

pressure,  P.  The  specific  expression  is 

ep+li;=D,(P*+T*)‘^^  (3) 


where  Di  and  D2  are  constants  and  P*  and  T*  are  as  defined  previously.  As  is  evident  from 
Equation  3,  the  concrete  material  cannot  undergo  any  plastic  strain  at  P*  =  -T*  without 
fracturing,  and  alternatively,  the  plastic  strain  to  fracture  increases  as  P*  increases.  A  third 
damage  constant,  (£p  +  Pp)„in ,  is  provided  to  allow  for  a  finite  amount  of  plastic  strain  to 

fracture  the  material.  This  is  included  to  suppress  fracture  from  low  magnitude  tensile  waves. 
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Figure  40.  Description  of  the  Holmquist-Johnson-Cook  (HJC)  Concrete  Model 
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Damage  due  to  plastic  volumetric  strain  is  included  in  Equations  2  and  3  because  the  concrete 
will  lose  cohesive  strength  during  air  void  collapse.  Under  most  circumstances,  the  majority  of 
the  damage  will  occur  from  equivalent  plastic  strain. 

The  hydrostatic  pressure-volume  relationship  is  presented  in  the  lower  right  comer  of  Figure  40. 
This  pressure-volume  response  is  separated  into  three  response  regions.  The  first  region  is  linear 
elastic  and  occurs  at  P  <  Pcmsh,  where  Pcrush  and  penish  are  the  pressure  and  volumetric  strain  that 
occur  in  a  uniaxial  stress  compression  test,  and  T  is  as  previously  defined.  The  elastic  bulk 

modulus  is  BCelastic  ~  Pcrush/ftcrush- 

The  second  region  is  referred  to  as  the  transition  region  and  occurs  at  Pcmsh  <  P  <  Piock-  In  this 
region,  the  air  voids  are  gradually  compressed  out  of  the  concrete  producing  plastic  volumetric 
strain.  Unloading  in  this  region  occurs  along  a  modified  path  that  is  interpolated  from  the 
adjacent  regions. 

The  third  region  defines  the  relationship  for  fully  dense  material  (all  air  voids  removed  from  the 
concrete).  The  air  voids  are  completely  removed  from  the  material  when  the  pressure  reaches 
Piock  and  the  relationship  is  expressed  as 

p  =  K,p:-i-K2p:^-HK3p:^  (4) 


where 


^  1  +  ^lock 


(5) 


The  modified  volumetric  strain,  p, ,  is  used  so  that  the  constants  (Ki,  K2,  and  K3)  are  equivalent 

to  those  used  for  material  with  no  voids.  The  standard  volumetric  strain  for  this  model  is 
p.  =  p  /  Po  —  1  for  current  density  p  and  initial  density  pQ.  The  locking  volumetric  strain  is 
M’lock  ~  Pgrain  ^  Po  ~  I  whcrc  pg,ai„  is  thc  grain  density.  This  is  identical  to  the  density  of  the 

material  with  no  air  voids. 

For  tensile  pressure,  P  =  in  the  elastic  region,  P  =  K,jl  in  the  fully  dense  region,  and  the 

pressures  are  interpolated  in  the  transition  region.  The  interpolation  factor  is 
F  =  (Pmax  -  Pcmsh)/(Ppiock  -  Pcmsh)  where  Pmax  is  the  maximum  volumetric  strain  reached  prior  to 
unloading  and  ppiock  is  the  volumetric  strain  at  Piock-  A  similar  method  is  used  for  compressive 
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unloading  except  that  the  higher  order  terms  andKjP^)  are  included.  The  tensile 

pressure  is  limited  to  T  (1  -  D). 

An  addition  to  the  model,  subsequent  to  Reference  58,  is  that  the  shear  modulus  varies  in  a 
manner  proportional  to  the  cmrent  bulk  modulus,  which  is  identical  to  using  a  constant  Poisson’s 
ratio. 

6.4  LIQUID  MATERIALS 

For  liquid  materials  the  three  normal  stress  components  are  expressed  as 


=  2(16, -(P  +  Q) 

(1) 

=  2m«,-(P  +  Q) 

(2) 

=  2ne.-(P  +  Q) 

(3) 

where  p  is  the  viscosity  coefficient  and  ,  Ey ,  and  are  the  normal  strain  rates.  P  is  the 

pressure  and  Q  is  the  artificial  viscosity  as  described  in  subsections  6.1.10  and  6. 1 . 1 1 .  For 
liquids  the  pressure  is  usually  limited  to  compression  only. 

The  three  shear  stress  are 


'Cxy=llYxy 

(4) 

II 

(5) 

'tyz=liYyz 

(6) 

where  ,  y^ ,  and  yy^  are  the  shear  strain  rates. 

6.5  BRITTLE  MATERIALS 

There  are  two  models  available  for  brittle/ceramic  materials,  the  original  Johnson-Holmquist 
(JH-1)  model  and  the  improved  Johnson-Holmquist  (JH-2)  model.  Both  models  were  developed 
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for  ceramic  materials,  but  they  can  also  be  used  for  other  brittle  materials  such  as  glass.  As  was 
the  case  for  the  crushable/concrete  models  imder  subsection  6.3,  the  strength,  pressure,  and 
damage  act  together  to  form  a  single  model. 

6.5.1  Johnson-Holmquist  JH’I  Model 

The  original  Johnson-Holmquist  brittle/ceramic  material  model  (Reference  59)  is  summarized  in 
Figure  41.  The  available  strength  (equivalent  stress),  a,  is  dependent  on  the  pressure,  P,  the 
dimensionless  strain  rate,  £*  =  £/£„  (for  £„  =  1.0  s"' ),  and  the  damage,  D.  For  undamaged 

material,  D  =  0;  for  partially  damaged  material,  0  <  D  <  1.0;  and  for  totally  damaged  (fractured) 
material,  D  =  1.0.  Note  that  the  strength  is  significantly  reduced  for  fractured  material  (D  =  1.0). 

T  is  the  maximum  tensile  hydrostatic  pressure  the  material  can  experience,  and  Si  and  S2  are  the 
strengths  of  the  intact  material  (for  £*  =  1.0)  at  the  compressive  pressures  Pi  and  P2,  respectively. 
After  the  material  has  fractured  (D  =  1.0),  the  slope  of  the  strength  is  given  by  Ce,  and  the 
maximum  fracture  strength  is  S3  (for  £*  =  l.O) . 

The  strain  rate  constant  is  C3.  If  Oo  is  the  available  strength  at  £*  =  1.0,  then  the  strength  at  the 
other  strain  rates  is 

a  =  a„(l  +  C3lnE*)  (1) 

It  can  be  seen  that  the  strength  increases  significantly  with  pressure,  which  is  consistent  with  the 
well-known  fact  that  brittle  materials  are  much  stronger  in  compression  than  they  are  in  tension. 
The  constants,  T,  Si,  and  Pi  can  generally  be  determined  from  quasi-static  and/or  dynamic 
(Hopkinson  bar)  tension,  torsion,  or  compression  tests;  and  the  strain  rate  constant,  C3,  can  be 
determined  from  comparable  quasi-static  and  dynamic  (Hopkinson  bar)  tests. 

The  higher  pressure  constants,  S2  and  P2,  generally  require  plate  impact  tests.  The  interpretation 
of  these  tests  can  be  difficult  because  generally  only  the  net  uniaxial  stress  can  be  measured.  To 
accurately  obtain  the  constants,  it  is  necessary  to  determine  both  the  hydrostatic  and  deviatoric 
components  of  stress. 


E_t19616.doc 


168 


STRENGTH 


e*>  1.0 
e*=  1.0 


The  post-fractxire  constant,  Ce,  can  be  bounded  by  two  different  types  of  tests.  A  lower  botmd 
can  be  established  by  axial  compression  testing  of  a  powdered  material  which  has  radial  pressure 
confinement,  and  an  upper  bound  can  be  established  by  axial  compression  testing  of  intact 
material  (using  displacement  control)  with  radial  pressure  confinement.  For  the  latter  technique. 
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it  is  the  strength  after  fracture  which  is  of  interest.  The  maximum  strength  of  the  fractured 
material,  S3,  may  sometimes  be  obtained  from  these  tests,  but  it  usually  requires  plate  impact 
testing. 

The  damage  for  fracture  is  accumulated  in  a  manner  similar  to  that  used  in  the  Johnson-Cook 
fracture  model,  presented  in  subsection  6.1.12.  It  is  expressed  as 

D  =  ZAep/e^  (2) 

where  ASp  is  the  plastic  strain  during  a  cycle  of  integration  and  =  f(P)  is  the  plastic  strain  to 

fracture  under  a  constant  pressure,  P.  Referring  to  Figure  41,  the  material  cannot  undergo  any 
plastic  strain  at  the  maximum  hydrostatic  tension,  T,  but  it  increases  to  £p  =  8^3,;  at  a 

compressive  pressure  of  P  =  DPi. 

The  hydrostatic  pressure  before  fracture  (D  <  1.0)  is  simply 

P  =  K,|i+K2ji"+K3p'  (3) 

where  Ki,  K2,  and  K3  are  constants  (Ki  is  the  bulk  modulus);  and  p  =  p/po  -  1  for  current  density 
p  and  initial  density  po.  For  tensile  pressures  (p  <  0),  Equation  3  is  replaced  by  P  =  KjP . 

Energy  effects  are  not  included. 

After  fracture  (D  =  1.0),  bulking  (pressure  increase  and/or  volumetric  strain  increase)  can  occur. 
Now  an  additional  incremental  pressure,  AP,  is  added,  such  that 

P  =  K,p  +  K2pVK3P^  +  AP  (4) 

The  pressure  increment  is  determined  from  energy  considerations.  Looking  back  to  the  strength 
model  in  Figure  41,  there  is  a  drop  in  strength  when  the  material  goes  from  an  intact  state  (D  < 
1.0)  to  a  fractured  state  (D  =  1.0).  This  represents  a  loss  in  the  elastic  internal  energy  of  the 
deviator  and  shear  stresses.  The  general  expression  for  this  internal  energy  is 

U  =  gV6G  (5) 

where  cr  is  the  equivalent  plastic  flow  stress  and  G  is  the  shear  modulus  of  elasticity. 
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The  loss  in  this  elastic  internal  energy  can  be  expressed  as 

AU  =  Ui-Uf  (6) 

where  Ui  is  the  elastic  energy  of  the  intact  material  before  fracture  (D  <  1 .0)  and  Uf  is  the  elastic 
energy  immediately  after  fi'acture  (D  =  1.0). 

This  energy  loss  (of  deviator  and  shear  stresses)  can  be  converted  to  potential  hydrostatic  internal 
energy  by  adding  AP.  An  approximate  equation  for  the  energy  conservation  is 

AP^f  +  AP  V  (2K, )  =  pAU  (7) 

where  pf  is  p  at  fracture  and  p  is  the  fraction  (0  <  p  <  1.0)  of  the  elastic  energy  loss  converted  to 
potential  hydrostatic  energy. 

The  first  term  (APp^)  is  the  approximate  potential  energy  for  p  >  0,  and  the  second  term 
[aP^  /  (2K,)]  is  the  corresponding  potential  energy  for  p  <  0. 

Solving  for  AP  gives 

AP  =  -K,Pf  +  ^(K,Pf)'  +2pK,AU  (8) 

Note  that  AP  =  0  for  p  =  0,  and  that  AP  increases  as  AU  increases  and/or  pf  decreases. 

Various  features  of  this  model  can  be  illustrated  by  the  three  examples  shown  in  Figure  42.  The 
tensile  strength  of  the  material  is  T  =  0.2  GPa,  the  imconfmed  compressive  strength  is  Si  =  S2  = 
2.0  GPa,  the  modulus  of  elasticity  is  E  =  220  GPa  and  Poisson’s  ratio  is  v  =  0.22.  Because  the 
height  is  H  =  1.0m  and  the  area  is  A  =  1.0m^  the  deflection  is  5  =  -Sz  and  the  force  is  F  =  -Oz. 

For  all  three  cases,  tiie  force,  F,  is  slowly  applied  until  6  =  0.02m,  and  then  it  is  slowly  released 
until  F  =  0.  The  paths  are  shown  for  strength  versus  pressure  and  for  force  versus  deflection. 
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Figure  42.  Examples  of  Material  Responses  with  the  JH-1  Ceramic  Model 
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For  Case  A,  the  material  cannot  develop  any  plastic  strain  =  o)  or  any  strength  after  fracture 

(of  =  0  for  C6  =  0).  All  of  the  elastic  energy  loss  (of  deviator  and  shear  stresses)  is  converted  to 
potential  hydrostatic  energy  (P  =  1.0).  Because  there  is  no  plastic  work,  the  external  work  must 
also  vanish.  This  is  clearly  shown  in  the  force  versus  deflection  relationship  for  Case  A.  The 
pressure  jump  at  fracture  (AP  =  0.56  GPa)  provides  for  this  conservation  of  energy. 

Case  B  is  similar  to  Case  A  except  that  the  material  is  allowed  to  undergo  a  small  amoimt  of 
plastic  strain  (ep  =  O.OOs)  prior  to  firacture.  Even  though  the  elastic  energy  loss  at  fi'acture,  AU, 

is  equal  to  that  of  Case  A,  the  pressure  jump  (AP  =  0.37  GPa)  is  less  because  pf  is  greater  at 
fi'acture. 

Case  C  is  similar  to  Case  B  except  that  there  is  strength  after  the  material  has  fi'actured  (of = 

0.5  P  for  Ce  =  0.5).  Here,  the  pressure  jump  (AP  =  0. 1 5  GPa)  is  reduced  firom  that  of  Case  B 
because  the  elastic  energy  loss  at  fracture,  AU,  is  less.  The  loading/imloading  path  is  very 
complex.  Of  special  interest  is  the  elastic  unloading  between  points  5  and  7.  Between  points  5 
and  6,  the  axial  deviator  stress,  Sz,  is  in  compression;  but  between  points  6  and  7,  the  same 
deviator  stress  is  in  tension.  The  net  stress,  az,  which  includes  the  hydrostatic  pressure,  remains 
in  compression.  At  point  7,  the  elastic  unloading  is  complete  and  the  material  flows  plastically 
between  points  7  and  8. 

6.5.2  Johnson-Holmquist  JH-2  Model 

The  improved  Johnson-Holmquist  brittle/ceramic  model  (Reference  60)  is  summarized  in 
Figure  43.  The  original  Johnson-Holmquist  (JH-1)  model  for  brittle  materials  includes  pressure- 
dependent  strength,  damage  and  fi'acture,  significant  strength  after  fi'acture,  bulking,  and  strain 
rate  effects,  as  presented  in  subsection  6.5.1 .  After  the  model  was  first  reported,  it  was 
implemented  into  several  codes  and  applied  to  a  variety  of  applications.  It  soon  became  apparent 
that  there  were  several  concerns. 

The  JH-1  model  does  not  allow  for  gradual  softening,  and  some  materials  clearly  show  a  gradual 
softening  during  flyer  plate  impact  tests.  Also,  results  for  some  applications  are  very  sensitive  to 
the  constants  used  in  the  model  and  there  is  not  a  straightforward  process  available  to  determine 
accurate  constants.  Finally,  the  jump  conditions  between  fractured  material  (damage  =  D  =  1.0) 
and  intact  material  (D  <  1.0)  caused  some  problems  for  Eulerian  codes,  whereby  the  material 
could  tend  to  heal  itself  after  fracture  had  occurred. 
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The  new  model  was  developed  to  address  these  concerns.  Reference  61  provides  a  description  of 
a  procedure  to  obtain  constants  for  glass,  as  well  as  some  plate  impact  and  penetration 
computations. 

A  general  overview  of  the  JH-2  model  is  shown  in  Figure  43.  Although  it  is  very  similar  to  the 
JH-1  model,  there  are  several  differences. 

•  The  material  begins  to  soften  when  the  damage  begins  to  accumulate  (D  >  0).  This 
allows  for  gradual  softening  of  the  material  under  increasing  plastic  strain.  (The  JH-1 
model  does  not  soften  imtil  D  =  1.0,  and  then  the  softening  occurs  instantaneously.) 

.  The  strength  and  pressure  are  normalized  by  the  strength  and  pressure  components  of 
the  Hugoniot  Elastic  Limit  (HEL),  which  allows  for  many  of  the  constants  to  be 
dimensionless.  This  can  be  very  helpful  when  comparing  different  materials,  and  when 
estimating  constants  for  materials  which  have  an  insufficient  database  to  determine 
constants. 

•  The  strength  and  damage  are  analytic  functions  of  the  pressure  and  other  variables. 

This  allows  for  parametric  variation  of  the  constants  in  a  more  systematic  manner.  (The 
JH-1  model  uses  multiple  linear  segments.) 

•  The  strength  generally  is  a  smoothly  varying  function  of  the  intact  strength,  fracture 
strength,  strain  rate,  and  damage.  It  is  well-suited  for  implementation  into  Eulerian 
codes. 

Returning  to  Figure  43,  the  normalized  equivalent  stress  is 

CT*=a;-D(a;-CT;),  (i) 

where  a*  is  the  normalized  intact  equivalent  stress,  a]  is  the  normalized  fracture  stress,  and  D  is 
the  damage  (0  <  D  <  1 .0). 

The  normalized  equivalent  stresses  (o’,  cr*,  )  have  the  general  form 


O*  —  c/CTheLj 


(2) 
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NORMALIZED  PRESSURE,  P*  =  P/Pu 


DAMAGE 
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Figure  43.  Description  of  the  Improved  Johnson-Holmquist  (JH-2)  Ceramic  Model 


F  t19616.doc 


where  cr  is  the  actual  equivalent  stress  and  cthel  is  the  equivalent  stress  at  the  HEL. 

The  normalized  intact  strength  is  given  by 

a;=A(P*+T*f(l  +  Clne*),  (3) 

and  the  normalized  fracture  strength  is  given  by 

a;=B(P*)"' (l+ClnE*).  (4) 

Note  that  the  normalized  fracture  strength  can  be  limited  by  Gj-  ^  Gf  (max) .  This  optional 
fracture  strength  parameter  is  included  to  provide  more  flexibility  in  defining  the  important 
fracture  strength.  It  also  allows  the  user  to  use  the  same  fracture  strength  as  used  for  the  JH-1 
model  (when  M  =  1.0). 

The  material  constants  are  A,  B,  C,  M,  N,  and  G*f  (max) .  The  normalized  pressure  is  P*  = 

P/Phel,  where  P  is  the  actual  pressure  and  Phel  is  the  pressure  at  the  HEL.  The  normalized 
maximum  tensile  hydrostatic  pressure  is  T*  =  T/Phft  ,  where  T  is  the  maximum  tensile 
hydrostatic  pressure  the  material  can  withstand.  The  dimensionless  strain  rate  is  e  =  £  /  , 

where  e  is  the  actual  strain  rate  and  £„  =  1.0  s~*  is  the  reference  strain  rate. 

The  damage  for  fracture  is  accumulated  in  a  manner  similar  to  that  used  in  the  JH-1  model  in 
subsection  6.5.1,  and  the  Johnson-Cook  fracture  model  in  subsection  6.1.12.  It  is  expressed  as 

D  =  ZAEp/eJ,  (5) 

where  As?  is  the  plastic  strain  during  a  cycle  of  integration  and  =  f  (P)  is  the  plastic  strain  to 

fracture  imder  a  constant  pressure,  P.  The  specific  expression  is 

£j=D,(P*+T*)°^  (6) 

where  Dj  and  Dz  are  constants  and  P*  and  T*  are  as  defined  previously  in  Equation  3.  Again,  the 
material  cannot  undergo  any  plastic  strain  at  P*  =  -T*,  but  £p  increases  as  P*  increases. 
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A  physical  explanation  of  damage  and  fracture  is  shown  in  Figure  44.  If  the  material  is  held 
under  a  constant  pressure  and  then  subjected  to  a  straining  deformation  at  a  constant  strain  rate, 
the  damage  begins  to  accumulate  when  the  material  begins  to  flow  plastically  (at  a  =  Oj).  The 
material  then  begins  to  soften  (relative  to  the  intact  strength).  This  softening  could  be  related  to 
the  material  going  from  a  larger  particle  size  to  a  smaller  particle  size  rmder  increased  plastic 
strain.  When  the  material  is  completely  damaged  (D  =  1.0),  the  strength  does  not  decrease  with 
increased  plastic  strain  (at  cr  =  Of). 


38.vsd  °  Cp  at  P  =  Pq 

T1 9616  Equivalent  Plastic  Strain,  e'’ 


Figure  44.  Strength,  Damage,  and  Fracture  Under  a  Constant  Pressure  and  Strain 
Rate  for  the  JH-2  Model 

Unfortunately,  it  is  not  generally  possible  to  perform  these  tests  at  sufficiently  high  pressures  of 
interest.  As  a  result,  the  damage  functions  and  fracture  strength  must  be  inferred  from  other  data. 

The  hydrostatic  pressure  before  fracture  begins  (D  =  0)  is  simply 

P  =  K,^  +  K2liVK3|i^  (7) 

where  Ki,  K2,  and  K3  are  constants  (Ki  is  the  bulk  modulus);  and  p  =  p/po  -  1  for  current  density 
p  and  initial  density  po.  For  tensile  pressure  (p  <  0),  Equation  7  is  replaced  by  P  =  K,p. .  Energy 

effects  are  assumed  to  be  insignificant. 
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After  damage  begins  to  accumulate  (D  >  0),  bulking  (pressure  increase  and/or  volumetric  strain 
increase)  can  occur.  Now  an  additional  incremental  pressure,  AP,  is  added,  such  that 

P  =  K,p  +  K2^'+K3^'+AP.  (8) 

The  pressure  increment  is  determined  fi-om  energy  considerations;  it  varies  from  AP  =  0  at  D  =  0 
to  AP  =  APmax  at  D  =  1 .0.  Figure  45  shows  how  AP  increases  as  D  increases.  The  incremental 
internal  elastic  energy  decrease  (due  to  decreased  shear  and  deviator  stresses)  is  converted  to 
potential  internal  energy  by  incrementally  increasing  AP.  The  decrease  in  the  shear  and  deviator 
stresses  occurs  because  the  strength,  a,  decreases  as  the  damage,  D,  increases,  as  shown  in 
Figure  43  and  Equation  1. 


P  =  +  KjpS  +  AP 
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Figure  45.  Description  of  incremental  Bulking  Pressure  for  the  JH-2  Modei 

The  general  expression  for  the  elastic  internal  energy  of  the  shear  and  deviator  stresses  is 

U  =  0^/60  (9) 

where  a  is  the  equivalent  plastic  flow  stress  and  G  is  the  shear  modulus  of  elasticity. 

The  incremental  energy  loss  is 
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AU  —  UD(t)  —  UD(t  +  At) 


(10) 


where  and  Ud(i+ at)  are  computed  from  Equation  9  using  CT,+At  for  both  energies. 

If  the  energy  loss,  AU,  is  converted  to  potential  hydrostatic  energy  through  AP,  an  approximate 
equation  for  this  energy  conservation  is 

(aP.,„  -AP.)p„^  +(aP.1,.  -AP.^)/2K,  =|3AU.  (11) 

The  first  term  [(aP,^.^,  -  AP,  ]  is  the  approximate  potential  energy  for  p  >  0  and  the  second 
term  [(aP,^^.^,  -  AP,^)/2K,]  is  the  corresponding  potential  energy  for  p  <  0. 

Solving  for  the  updated  AP  gives 

AP,..  +i/(K,n,..+AP,f+2|3K,AU.  (12) 

As  was  the  case  for  the  JH-1  model,  the  JH-2  model  also  gives  AP  =  0  when  P  =  0,  where  p  is  the 
fraction  (0  <  p  <  1)  of  the  elastic  energy  loss  converted  to  potential  hydrostatic  energy. 

Various  features  of  the  JH-2  model  can  be  illustrated  by  the  three  examples  shown  in  Figure  46. 
The  tensile  strength  of  the  material  is  T  =  0.2  GPa,  the  HEL  =  2.79  GPa,  the  intact  strength  is 
assumed  to  be  a*  =  0.93  (P*  -i-  T*)°  * ,  the  modulus  of  elasticity  is  E  =  220  GPa,  and  Poisson’s 

ratio  is  V  =  0.22.  There  is  no  strain  rate  effect  so  that  C  =  0  in  Equations  3  and  4.  From  the 
preceding  constants,  the  following  additional  constants  can  be  obtained: 


'''  =  3(l-2v)  = 

(13) 

(14) 

Phf.t.  “  KiM'hel  ^2M'hel  ^3^HEL  ”  GPa , 

(15) 

^HEL  — 

(16) 
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Figure  46.  Examples  of  Material  Responses  with  the  JH>2  Ceramic  Model 
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where  HEL  is  the  net  axial  stress  for  the  Hugoniot  Elastic  Limit,  Phel  is  the  pressure  component 
of  the  HEL,  and  Phel  =  0.01 1 17  is  the  corresponding  p  at  Phel  (for  K2  =  K3  =  0).  For  a  known 
HEL  and  G,  Phel  can  be  obtained  iteratively  from 


HEL  —  K.jPjj£l  K2M'hel  ^sM'hel  3 


l^HEL 


(17) 


Because  the  height  is  H  =  1  .Om  and  the  area  is  A  =  1  .Om^,  the  deflection  is  5  =  -8z  and  the  force 
is  F  =  For  all  three  cases,  the  force,  F,  is  slowly  applied  imtil  6  =  0.05m,  and  then  it  is 
slowly  released  until  F  =  0.  The  paths  are  shown  for  strength  versus  pressure  and  for  force 
versus  deflection. 

For  Case  A,  the  material  cannot  develop  any  plastic  strain  (sp  =  oj  or  any  strength  after  fracture 

(Of  =  0) .  All  of  the  elastic  energy  loss  (of  deviator  and  shear  stresses)  is  converted  to  potential 

hydrostatic  energy  (P  =  1.0).  Because  there  is  no  plastic  work,  the  external  work  must  also 
vanish.  This  is  shown  in  tide  force  versus  deflection  relationship  for  Case  A.  The  pressure  jump 
(AP  =  0.56  GPa)  provides  for  this  conservation  of  energy. 

Case  B  is  similar  to  Case  A  except  the  material  is  allowed  to  accumulate  some  plastic  strain 
during  fracture.  Taking  Dj  =  0.005  and  D2  =  1.0  in  Equation  6  gives  a  pressure  dependent 
fracture  strain,  .  It  can  be  seen  in  Figure  46  that  the  material  gradually  softens  between  point  2 

(where  the  damage  begins  to  accumulate)  and  point  3  (where  the  damage  is  complete  and  there  is 
no  strength).  The  bulking  pressure,  AP,  is  also  generated  between  these  two  points. 

Case  C  is  similar  to  Case  B  except  that  there  is  strength  in  the  fractured  material 

[of  =  0.31  (P*)”  *] .  The  constants  from  Equation  4  are  assumed  to  be  B  =  0.31,  m  =  0.6,  and  C  = 

0.  Here  the  response  is  very  complex,  with  the  primary  regions  of  interest  as  follows: 

From  Points  1  to  2  the  material  loads  elastically  rmtil  the  intact  strength,  ai,  is  encoimtered  at 
point  2.  From  Points  2  to  3  the  material  flows  plastically,  moving  from  the  intact  strength  at 
point  2  to  the  fracture  strength  at  point  3.  The  damage  goes  from  D  =  0  at  point  2  to  D  =  1 .0  at 
point  3.  Similarly,  the  bulking  pressure  goes  from  AP  =  0  at  point  2  to  AP  =  0.65  GPa  at  point  3. 
From  Points  3  to  4  the  material  continues  to  flow  plastically  along  the  fracture  strength  envelope, 
af.  From  Points  4  to  5  the  loading  direction  is  reversed  at  point  4  and  the  material  rmloads 
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elastically.  At  point  4,  the  axial  deviator  stress  is  in  compression.  At  point  5,  the  axial  deviator 
stress  is  completely  unloaded  such  that  the  resulting  stress  is  due  only  to  the  hydrostatic  pressure. 
From  Points  5  to  6  the  elastic  unloading  continues  and  the  axial  deviator  stress  goes  into  tension. 
The  net  axial  stress  remains  in  compression,  however,  because  the  pressiue  is  greater  than  the 
deviator  stress.  The  fracture  stress,  Gf,  is  encoimtered  at  point  6.  From  Points  6  to  7  the  material 
unloads  plastically  along  the  fracture  strength  envelope.  At  point  7,  both  the  pressure  and  the 
axial  deviator  stress  (strength)  go  to  zero. 

6.6  REACTIVE  EXPLOSIVE  MATERIALS 

For  the  Reactive  Explosive  materials  the  bum  fraction  is  determined  as  a  function  of  the  pressure 
history.  For  the  Explosive  Materials  (described  previously  in  subsection  6.2)  the  bum  fraction 
could  be  predetermined  based  on  the  detonation  velocity  and  the  point  of  detonation.  The 
Reactive  Explosive  models  generally  require  a  very  fine  grid  and  significantly  more  CPU  time. 

6.6.1  Tarver  Ignition  and  Growth  Model 

The  Ignition  and  Growth  model  is  presented  in  Reference  62.  This  reference  uses  the  JWL 
Equation  of  State  (presented  previously  in  subsection  6.2.2)  for  both  the  pressure  in  the  solid 
material  and  the  gas  products.  The  bum  rate  fraction,  F,  goes  from  F  =  0  (solid  only,  no  reaction) 
to  F  =  1.0  (complete  reaction).  The  bum  rate  is 

F  =  aF/at  =  l(l-F)®(p/po-l-A)’‘+G,(l-F)''F'’P''+G2(l-F)®F°P^  (1) 

The  three  basic  terms  correspond  to  ignition,  growth,  and  completion.  F  is  the  mass  fraction  of 
the  explosive  which  has  reacted,  t  is  the  time,  p  is  the  current  density,  po  is  the  initial  density,  and 
P  is  the  pressiire.  The  constants  are  I,  Gi,  G2,  A,  B,  C,  D,  E,  G,  X,  Y,  and  Z.  The  dimensionless 
constant.  A,  prohibits  ignition  until  (p/po  -  1  -  A)  >  0. 

The  bum  rate,  F,  is  given  by  Equation  1,  and  the  bum  firaction,  F,  is  numerically  integrated  as 
follows: 

pt+At  ^F‘+F'^At  (2) 
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where  At  is  the  integration  time  increment  and  F‘^  is  based  on  the  volumetric  strain  at  time  =  t  + 
At  and  the  pressure  at  time  =  t.  F  =  0  indicates  all  solid  behavior  (no  detonation),  and  F  =  1 .0 
indicates  all  gas  products  behavior  (complete  detonation).  The  bum  fraction  is  an  element 
variable  which  must  be  carried  from  cycle  to  cycle. 

There  are  two  other  element  variables  which  must  also  be  carried  from  cycle  to  cycle.  VFg  is  the 
fraction  of  the  element  volume  which  contains  gas  products,  and  ETg  is  the  total  internal  energy 
which  has  been  added  to  the  gas  products. 


ETg  =  (Eg-E„)V„F 


(3) 


where  Eg  is  the  internal  energy  per  initial  volume  in  the  gas  products,  Eo  is  the  initial  internal 
energy  per  initial  voliune  in  the  explosive,  Vo  is  the  initial  volume  of  the  element,  and  F  is  the 
bium  fraction. 

The  initial  estimate  for  VFg  is  given  by 

VF‘^^‘ =VF‘+[aF/(1-F')](i-VF^)  (4) 

where  AF  =  F*^^*  —  F‘  is  the  incremental  bum  fraction  between  cycles.  The  quantity  AF/(1  —  F') 
is  simply  the  fraction  of  unbumed  mass  which  was  burned  during  the  past  cycle.  Similarly,  the 
quantity  (l  -  VFg  )  is  the  volume  fraction  of  imbumed  material.  This  equation  essentially 

converts  the  volume  of  solid  burned  during  the  past  cycle  to  an  equivalent  voliune  of  gas 
products.  Equation  4  is  only  an  estimate  and  is  therefore  altered  in  subsequent  iterations. 

The  relative  volumes  for  the  entire  element,  solid  portion,  and  gas  products  portion,  are 
designated  as  V,  Vg,  and  Vg ,  respectively.  The  all  have  the  form  V  =  V  /  V^,  where  V  is  the 

current  volume  and  Vo  is  the  initial  volume. 

The  relative  volumes  of  the  solid  and  gas  products  portions  at  time  =  t,  are  given  by 

V;  =  V*  (l  -  VFg*  -  AVf)  /  (l  -  F*^"" )  (5) 

Vg*  =  V*  (vFg*  +  AVf)  /  F*^^*  (6) 
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where  AVF  =  [aF  /  (l  -  F*  )]^1  -  VFg )  is  required  to  represent  the  volume  of  material  which 
burned  during  the  previous  cycle. 

Similarly,  the  relative  voliimes  at  time  =  t  +  At,  are 


•wt+At  _  -^t+At 
Vg  V 

(l-VFg‘^^')/(l-F*'"'^*) 

(7) 

“^t+At  _ ^t+At 

ypt+At  !  pt+At 

(8) 

For  the  first  iteration,  VFg^^'  is  given  by  Equation  4;  but  in  subsequent  iterations,  it  is  varied. 
The  corresponding  volumetric  strain  rates  are 


(9) 

e_=(7'“-V.‘)/At 

(10) 

The  internal  energy  distributions  must  also  be  determined  before  the  JWL  Equation  of  State 
(EOS)  can  be  applied. 

The  total  added  internal  energy  is  given  by 

et=(e-eJv„  (11) 

where  E  is  the  cmrent  internal  energy  per  initial  volume  and  Eo  is  the  initial  internal  energy  per 
initial  volume  in  the  explosive. 

The  total  added  energy  in  the  gas  products,  ETg,  is  given  in  Equation  3.  The  total  added  energy 
in  the  solid  is  then 


ET,  =  ET-ETg  (12) 

As  the  explosive  bums,  the  internal  energy  of  the  mass,  which  changes  fi-om  solid  to  gas  products 
during  an  integration  cycle,  must  be  transferred  from  the  solid  to  the  gas  products. 
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ET,'^  =ET3‘-AET 


(13) 


ETg*^  =ET^+AET 


(14) 


where  the  energy  change  is  given  by 
AET  =  ET3AF/(1-F‘) 


(15) 


Now  the  internal  energies,  per  initial  volume,  for  the  JWL  EOS  can  be  determined. 

E.=E‘+ET./[V,-  (l-F)]+AEj  (16) 

E,  =E„  +  ET,/(v,f)  (17) 

where  E®  and  E^  are  the  initial  internal  energies  per  initial  volume,  for  the  solid  and  the  gas 
products.  In  Equation  16,  AEd  is  the  incremental  internal  energy,  generated  during  the  past  cycle, 
due  to  the  material  strength  (shear  and  deviator  stresses)  of  the  solid  material.  (As  described  in 
subsection  6.1.9). 


The  JWL  EOS  is  given  in  subsection  6.2.2.  After  the  pressures  for  both  the  solid  and  gas 
products  have  been  computed,  an  equilibrium  pressure  difference  is  computed  by 


AP  = 


ip.-pj 

ip.i  +  IP, I 


(18) 


If  AP  <  AP3„o^^  ,  then  equilibriiun  has  been  attained.  For  AP  >  AP^jio^^j ,  a  new  value  of  VFg  is 

assumed  (using  a  binary  iteration)  and  the  previous  process  is  repeated.  Limited  experience 
indicates  reasonable  results  are  obtained  using  APguo^gj  =  0.01 . 


Although  Reference  62  (and  the  discussion  in  this  subsection)  uses  the  JWL  EOS  for  both  the 
solid  material  and  the  gas  products,  it  is  possible  to  use  other  Equations  of  State  with  the  bum 
rate  given  in  Equation  1. 
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6.6.2  Mader  Forest  Fire  Model 

The  Forest  Fire  model  is  presented  in  Reference  63.  Here  the  bum  rate  is  expressed  as 

^  =  (l-F)exp(Co+C,P+C,P^+ ...  +C„P")  (1) 

where  Cq.  .  .Cm  are  material  constants  and  P  is  the  pressure.  The  bum  fraction,  F,  is  as  defined  in 
subsection  6.6.1.  F  =  0  indicates  all  solid  behavior  (no  detonation),  and  F  =  1.0  indicates  all  gas 
products  behavior  (complete  detonation).  The  numerical  implementation  is  similar  to  that 
described  in  subsection  6.6.1. 

6.7  RDG  STRENGTH  AND  FRACTURE  MODEL 

The  RDG  strength  and  fracture  model  is  described  in  References  64  and  65.  This  model  uses 
either  the  Johnson-Cook  or  Bodner-Partom  strength  model  and  the  Mie-Gruneisen  Equation  of 
State.  These  are  coupled  together  with  a  void  nucleation  fracture  model.  The  current  strength 
and  pressure  are  determined  as  a  function  of  the  intact  material  properties  and  the  void  volume 
fraction.  The  specific  algorithm  is  not  presented  in  this  report.  This  is  a  complex  model,  and  the 
References  should  be  consulted  for  a  detailed  description. 
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